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This paper is a review of results which have been recently obtained by applying mathematical 
concepts drawn, in particular, from differential geometry and topology, to the physics of Hamiltonian 
dynamical systems with many degrees of freedom of interest for statistical mechanics. 

The first part of the paper concerns the applications of methods used in classical differential 
geometry to study the chaotic dynamics of Hamiltonian systems. Starting from the identity between 
the trajectories of a dynamical system and the geodesies in its configuration space, when equipped 
with a suitable metric, a geometric theory of chaotic dynamics can be developed, which sheds new 
hght on the origin of chaos in Hamiltonian systems. In fact, it appears that chaos can be induced not 
only by negative curvatures, as was originally surmised, but also by positive curvatures, provided 
the curvatures are fluctuating along the geodesies. In the case of a system with a large number of 
degrees of freedom it is possible to approximate the chaotic instability behavior of the dynamics by 
means of a geometric model independent of the dynamics, which allows then an analytical estimate 
of the largest Lyapunov exponent in terms of the averages and fluctuations of the curvature of the 
conflguration space of the system. 

In the second part of the paper the phenomenon of phase transitions is addressed and it is here that 
topology comes into play. In fact, when a system undergoes a phase transition, the fluctuations of 
the conflguration-space curvature, when plotted as a function of either the temperature or the energy 
of the system, exhibit a singular behavior at the phase transition point, which can be qualitatively 
reproduced using geometric models. In these models the origin of the singular behavior of the 
curvature fluctuations appears to be caused by a topological transition in configuration space, which 
corresponds to the phase transition of the physical system. This leads us to put forward a Topological 
Hypothesis (TH). The content of the TH is that phase transitions would be related at a deeper level 
to a change in the topology of the conflguration space of the system. We will illustrate this on a 
simple model, the mean-fleld XY model, where the TH can be checked directly and analytically. 
Since this model is of a rather special nature, namely a mean-fleld model with infinitely-ranged 
interactions, we discuss other more realistic (non mean-field-like) models, which can not be solved 
analytically, but which do supply direct supporting evidence for the TH via numerical simulations. 
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I. INTRODUCTION 

This paper deals with the application of concepts drawn from mathematics, in particular from differential geometry 
and topology, to problems in statistical physics. The mathematical tools involved come from Riemannian geometry 
and from Morse theory, respectively. As to the physics, the applications of these concepts will be brought to bear on 
dynamical systems with many degrees of freedom, including eventually the thermodynamic limit. 

In order to contain this report to a reasonable size and yet make it accessible to as wide a readership as possible, 
and since it makes use of concepts which might not be known to everyone, we chose the following format. 
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The first part of the main text is aimed at a reader who is famihar with the basics of Riemannian geometry, for 
example at the level of a course in general relativity. As to the second part, the knowledge of Morse theory at an 
elementary level is assumed. However, for those who are not familiar with these branches of mathematics, we have 
provided in extended Appendices the main points which are needed to follow the exposition. Similarly, we assume that 
the reader is familiar with the basics of dynamical systems theory, but again we summarized in an Appendix the main 
concepts. In all cases references to the literature are made for the details. Both the main text and the Appendices 
are written as a compromise between mathematical rigour and a physcist's accessibility; in case of conflict, we always 
favoured the latter. 

This way wc hope that a reader familiar with the basic mathematical tools will be able to read the report straight- 
forwardly. We have made a special effort to emphasize logical structure and physical content, and we hope that the 
report will provide a clear survey of what has been achieved applying geometrical methods to dynamical systems and 
statistical mechanics. At the same time wc would like this paper to allow the reader to familiarize herself or himself 
with this new field and to stimulate new developments and contributions to the many points which are still open. 

Throughout the paper wc will consider classical Hamiltonian dynamical systems with N degrees of freedom, confined 
in a finite volume^, whose Hamiltonian is of the form 



where the q's and the p's are, respectively, the coordinates and the conjugate momenta of the system. Our emphasis is 
on systems with a large number of degrees of freedom. The dynamics of the system (1) is defined in the 2iV-dimensional 
phase space spanned by the q^s and the p's. Our aim is to relate the dynamical and the statistical properties of the 
system (1) with the geometrical and topological properties of the phase space where the dynamical trajectories of the 
system live. It turns out that as long as we consider Hamiltonians of the form (1) we can restrict ourselves to the 
study of the geometry and the topology of the A/'-dimensional configuration space^ without loosing information. In 
fact as we shall sec in Sec. HI the dynamical trajectoricis can be seen as geodesies of the configuration space, provided 
the latter has been endowed with a suitable metric. As to the topology, later on we shall see that also all the relevant 
information on the topology of phase space, from the point of view of Morse theory, is encoded in the potential energy 
function V{qi, . . . , (/at), so that also the topological investigation can be restricted to the configuration space. 

This is similar to what happens in the classical statistical mechanics of Hamiltonian systems of the form (1), where 
the momenta can be integrated out and the statistical measure can be defined on the configuration space alone. We 
remark that this is true for both the microcanonical and the canonical ensemble. 

The structure of the paper is the following: after a short historical introduction (Sec. II), the main body of the 
paper is organized in two parts, which, though being tightly related, are also to a great extent independent from each 
other, so they could be also read separately without encountering too many difficulties. 

The first part (Sees. Ill and IV) concerns the applications of tools belonging to classical differential geometry to 
study the chaotic dynamics of Hamiltonian systems of the form (1), in particular for those with a large number of 
degrees of freedom. Starting from the identity between the trajectories of the dynamical system (1) and the geodesies 
of the configuration space equipped with the Jacobi or the Eisenhart metric, we develop a geometric theory of chaotic 
dynamics which sheds a new light on the origin of chaos in Hamiltonian systems. In fact it turns out that chaos can 
be induced not only by negative curvatures, but also by positive curvatures, provided the curvatures are fiuctuating 
along the geodesies. In the case of a large number of degrees of freedom it is possible to describe the instability of the 
dynamics by means of a geometric model independent of the dynamics, which provides an analytical estimate of the 
largest Lyapunov exponent in terms of the averages and fluctuations of the curvature of the configuration space. We 
clarified the exposition with the application of the general concepts to two special examples; however, since some of 
the calculations are rather lengthy, we did not provide all the details, referring the reader for those to the appropriate 
literature. 

While in the first part we deal with the application of differential geometry to dynamical systems with a large 
number of degrees of freedom, we do not touch upon one of the most spectacular properties of large systems, namely 
that in the thermodynamic limit N ^ oo they may show sharp phase transitions. In the second part (Sees. V 
and VI) we address this point, and it is here that topology comes into play. In fact, when a system undergoes 



'^We will mainly apply our results to systems defined on a lattice, so that we will not explicitly refer to the volume of the 
system. 

^Actually we will also consider an enlarged configuration space with two extra dimensions. 




(1) 
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a phase transition, the fluctuations of the configuration-space curvature, when plotted as a function of either the 
temperature or the energy, have a singular behavior at the transition point which can be qualitatively reproduced 
using a geometric model. In such a model the origin of the singular behavior of the curvature fluctuations resides 
in a topological change. This leads us to put forward a Topological Hypothesis (TH). The content of the TH is that 
phase transitions (at least, continuous phase transitions) would at a deeper level be related to a particular change in 
the topology of the conflguration space of the system. We will illustrate this on a simple model, the mean-field XY 
model, where the TH can be checked directly and analytically. Since this model is of a rather special nature, namely 
a mean-fleld model with infinitely-ranged interactions, we discuss other more realistic (non mean-field-like) models, 
which can not be solved analytically, but which do supply indirect supporting evidence for the TH via numerical 
simulations. 

We emphasize that the purpose of the work discussed in the second part of this Report is not to extend the existing 

theory of phase transitions for systems described by a classical Hamiltonian. Rather, we try to extend the foundation 
for the occurrence of phase transitions to transitions in the topology of the configuration space of the system undergoing 
a phase transition. This way we hope to make a new connection between a branch of pure mathematics (topology) 
and a branch of statistical mechanics (phase transitions). Such a connection appears to lead to a new approach to 
phase transitions. 

A brief historical summary, to place the content of the present review in its historical context, is given in the next 
Section. 



II. HISTORICAL REMARKS 

Without attempting to be exhaustive, a few historical comments might be helpful to place the recent contributions 
about the geometrical approach to dynamics and statistical physics which are reviewed in the present article, in a 
more general context. This makes the present Section an exception to the self-containedness of this review paper, 
because we mention here concepts which arc not necessary to understand the topics treated in the rest of the paper. 

The idea of looking at the collection of solutions of Newton's equations of motion from a geometric point of view 
dates back to H. Poincare and to the development of the qualitative theory of differential equations. Tackling the 
famous problem of the integrability of the three-body problem, Poincare also discovered that generic classical Hamil- 
tonian systems, in spite of their deterministic nature, lack predictability, i.e. were unstable, because of their extreme 
sensitivity to the initial conditions. Such an instability of classical dynamics originates in homoclinic intersections, 
which Poincare described in his Methodes Nouvelles de la Mecanique Celeste [1] without "even attempting to draw" 
them (see Appendix C). However his geometric treatment of dynamics, later developed by E. Cartan among others, 
involves submanifolds of phase space using what is now called symplectic geometry [2]. Although of undeniable ele- 
gance, symplectic geometry is not very helpful to advance our knowledge about those regions in phase space where 
the dynamics is unstable, knowledge which is relevant for statistical mechanics. 

The name of Poincare, together with that of E. Fermi, is also associated with an important theorem about the 
non-existence of analytic integrals of motion, besides energy, for generic nonlinear Hamiltonian systems describing at 
least three interacting bodies [1,3]; this is the origin of the accessibility of the whole constant energy hypersurface of 
phase space which is determined by the initial conditions. The Poincare-Fcrmi theorem has been generally considered 
by the physics community as sufficient to legitimate classical statistical mechanics from a dynamical viewpoint"*. 

One had to wait until the 1940s when a qualitatively new attempt c;nierged to make use of geometric conccipts in the 
investigation of newtonian dynamics and its connection with statistical mechanics. Though very different from the 
clear mathematical expositions of Poincare, and much more in the spirit of a physicist, it was N. S. Krylov [5] in his 
doctoral dissertation, who showed for the first time the existence of a close relationship between dynamical instability 
(seen as the exponential amplification of small deviations in the initial conditions of a collection of colliding objects 
representing idealized atoms in a gas) and phase space mixing. Phase mixing is a stronger property than ergodicity 
and is far more relevant to physics than ergodicity. In fact, while ergodicity assures the equality of time and phase 
space averages of physical quantities, phase mixing addresses the rate of approach of time to ensemble averages. Phase 
mixing entails thus the convergence of time averages to ensemble averages in a finite time. In modern terms, Krylov 



^Kormogorov-Arnord-Moscr(KAM) theory [4] might seem capable of explaining how, in spite of this "no-go" theorem due to 
Poincare and Fermi, ergodicity could fail. However, the exceedingly tiny - and lastly vanishing with the number A'^ of degrees 
of freedom - nonintegrable perturbations that arc required to keep a positive measure of regular regions in phase space, do, in 
general, not have any physically appreciable effect even when only a few bodies are considered. 
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realized the necessity of chaotic dynamics to obtain phase mixing and to make the connections between dynamics and 
statistical mechanics stronger. Moreover Krylov, in view of what we are going to discuss in the next Sections, also 
has the great historical merit of having attempted for the first time to bridge the dynamical foundations of statistical 
mechanics with a widely developed and powerful field of mathematics: Riemannian diflferential geometry. Krylov knew 
certain mathematical results concerning the properties of geodesic flows on compact negatively curved manifolds by 
Hadamard, Hedlund and Hopf [6], and he realized their potential interest to physics, once newtonian dynamics is 
rephrased in terms of Riemannian geometric language. Such a possibility was well known since the beginning of the 
century, mainly due to the work of T. Levi-Civita, in particular that the principle of stationary action entails the 
identity of a classical mechanical flow with a geodesic fiow in a configuration space, endowed with a suitable metric. 
Krylov's efforts concentrated on the analysis of the properties of physical systems which move in negatively curved 
regions in configuration space. For example, he discussed how the presence of an inflection point in the Lennard- Jones 
potential could influence the dynamics of a dilute gas (through the appearance of regions of negative scalar curvature 
in configuration space) and lead to a strong instability of the dynamics. These attempts have been very influential on 
the development of the so-called abstract ergodic theory [7], where Anosov flows [8] (e.g., geodesic flows on compact 
manifolds with negative curvature) play a prominent role. Ergodicity and mixing of these flows have been thoroughly 
investigated. To give an example, Sinai proved ergodicity and mixing for two hard spheres by just showing that such 
a system is equivalent to a geodesic flow on a negatively-curved compact manifold [9] 

From time to time, Krylov's intuitions have been worked out further by several physicists for whom we refer to Refs. 
[10-20]. They invariably discovered, much to their surprise, that geodesic flows associated with physical Hamiltonians 
do not live on negatively curved manifolds, despite their chaoticity even if the latter is well developed; only very few 
exceptions to this are known, in fact the two low dimensional models discussed in Refs. [12,15], where chaos is actually 
associated with negative curvature. Worse, for certain models the regions of negative curvature of the mechanical 
manifolds apparently shrink by increasing the number N of degrees of freedom. 

This somewhat biased search for negative curvatures has been the main obstacle to an effective use of the geometric 
framework originated by Krylov to explain the source of chaos in Hamiltonian systems. On the other hand, it is true 
that the Jacobi equation, which describes the (in)stability of a geodesic flow, is in practice only tractable on negatively 
curved manifolds, formidable mathematical difficulties are encountered in treating the (in)stability of geodesic flows 
on manifolds of non-constant and not everywhere negative curvature. Moreover, for this kind of problems, intuition 
can hardly help. However, the advent of computers has been here of invaluable help. In this connection, it may not 
be out of place to quote here some of the sentences which S. Ulam remembered from the far- looking conversations he 
had with E. Fermi and J. von Neumann [21]: 

"After the war, during one of his frequent smnmer visits to Los Alamos, Fermi became interested in the 
development and potentialities of the electronic computing machines. He held many discussions with me 
on the kind of future problems which could be studied through the use of such machines. We decided to 
try a selection of problems for heuristic work where in the absence of closed analytic solutions experimental 
work on a computing machine would perhaps contribute to the understanding of properties of solutions. 

(...) Fermi expressed often the belief that future fundamental theories in physics may involve nonlinear 
operators and equations, and that it would be useful to attempt practice in the mathematics needed for 
the understanding of nonlinear systems (...)." 

As a matter of fact, only during the last few years an interplay between analytic methods and numerical simulations 
has made possible to overcome the mentioned difficulties, proving the effectiveness of the Riemannian geometric 
approach to dynamical systems of interest to statistical mechanics, field theory, and condensed matter physics [22-39]. 
As we shall see in the following, this has extended the domain of application of geometric techniques, and has also 
introduced a new point of view about the origin of chaos in Hamiltonian systems, as well as new methods to describe 
and understand it. 

The use of Finsler manifolds (generalizations of Riemannian manifolds that allow the geometrization of velocity 
dependent potentials as well) has also been proposed in Refs. [40,41]. An analysis of dynamics based on the geometry 
of trajectories, rather than of the manifolds on which they move, has been proposed in Ref. [42] . 

For what concerns the use of geometric and topologic concepts in statistical mechanics, we must distinguish be- 
tween macroscopic and microscopic phase spaces. A macroscopic phase space is a low-dimensional space spanned 
by macroscopic variables, like temperature, pressure, volume, chemical potential etc, in other words it is in general 
a parameter space. In the '70s some applications to the study of phase transitions of the theory of singularities of 
differentiable maps (popularly known as Catastrophe theory), which includes Morse theory, were proposed. These 
followed R. Thom's remark that the simplest example of the classical critical point as it appears in the van der Waals 
equation corresponds to the Riemann-Hugoniot catastrophe [43,44]. 

An elegant and deep formulation of phase transitions related to structural instability and using one of the most 
beautiful theorems in difii'erential topology, the Atiyah-Singer index theorem, was proposed by M. Rasetti in Ref. [45]. 
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Other very recent proposals of geometric and topologic methods in macroscopic phase spaces have been put forward 
in Refs. [46,47]. 

In recent papers some elements of the geometry of constant energy hypersurfaces T.e oi phase space have been used 
for the microcanonical definitions of the temperature and the specific heat [48,49], in Ref. [50] a topological property 
of the E E has been related to their mean curvature from which a relationship between thermodynamics and topology 
emerged. 

The description of phase transitions through geometric and topologic changes in the microscopic phase space 
has never been considered until very recently. It appeared as a natural consequence of the above mentioned new 

developments in the Ricmannian theory of Hamiltonian chaos. These newly proposed ideas, as well as the conceptual 
path that led from the geometry of dynamics to topology and phase transitions, are reviewed in this paper. 

III. RIEMANNIAN GEOMETRY AND HAMILTONIAN DYNAMICS 

A Hamiltonian system whose kinetic energy is a quadratic form in the velocities is referred to as a natural Hamilto- 
nian system. Every Newtonian system, that is a system of particles interacting through forces derived from a potential, 
i.e. of the form (1), belongs to this class. The trajectories of a natural system can be seen as geodesies of a suitable 
Riemannian manifold. This classical result is based on a variational formulation of dynamics. In fact Hamilton's 
principle states that the natural motions of a Hamiltonian system are the extrema of the functional (Hamiltonian 
action <S) 




where C is the Lagrangian function of the system, and the geodesies of a Riemannian manifold are the extrema of the 
length functional 




where s is the arc-length parameter. Once a connection between length and action is established, by means of a 
suitable choice of the metric, it will be possible to identify the geodesies with the physical trajectories. 

A. Geometric formulation of the dynamics 

The Riemannian formulation of classical dynamics is far from unique, even if we restrict oureselves to the case of 
natural systems. There are many possible choices for the ambient space and its metric. The most commonly known 
choice — dating back to the nineteenth century — is the so-called Jacobi metric on the configuration space of the 
system. Actually this was the geometric framework of Krylov's work. Among other possibilities, we will also consider 
a metric originally introduced by Eiscnhart on an enlarged configuration space-time. The choice of the metric to be 
used will be dictated mainly by convenience. 

These choices certainly do not contain all the possibilities of geometrizing conservative dynamics. For instance, 
with regards to systems whose kinetic energy is not quadratic in the velocities — the classical example is a particle 
subject to conservative as well as velocity-dependent forces, such as the Lorentz force — it is impossible to give a 
Riemannian gcometrization, but becomes possible in the more general framework of a Finsler geometry [51]. However, 
we will not consider this here, and restrict ourselves to standard Hamiltonian systems. 

For a summary of the notation and the concepts of differential geometry that will be used in the following we refer 
the reader to Appendix A. The summation convention over repeated indices will be always used throghout the paper, 
if not explicity stated otherwise. 

1. The Jacobi metric 

Let us consider an autonomous dynamical system, i.e., a system with interactions which do not explicitly depend 
on time, whose Lagrangian can be written as 

jC=T-V= ^aijq'qi - V{q) , (4) 
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where the dot stands for a derivative with respect to the parameter on which the q's depend"*, and q is a. shorthand 
notation for all the coordinates qi, . . . , q^- Both these conventions will be used throughout the paper, when there is 
no possibility of confusion. 

The Hamiltonian Ti = T -\- V is &n integral of motion, whose value, the energy E, is a conserved quantity. Hence 
Hamilton's priciple can be cast in Maupertuis' form [2]: the natural motions of the system are the stationary paths 
in the configuration space M for the functional 

among all the isoencrgctic curves, i.e. the curves ^{t) connecting the initial and final points parametrized so that the 
Hamiltonian q) is a constant equal to the energy E. The fact that the curves must be isoenergetic with energy 
E implies that the accessible part of the configuration space is not the whole M, but only the subspace Me C M 
defined by 

ME = {qeM: V{q) < E} . (6) 

In fact a curve 7' which lies outside Me will never be parametrizable in such a way that the energy is E, because 7' 

will then pass through points where V > E and the kinetic energy is positive"''. 

The kinetic energy T is a homogeneous function of degree two in the velocities, hence Euler's theorem implies that 

^T = q^§, (7) 

and Maupertuis' principle reads as 

dA = S [ 2Tdt = . (8) 



The configuration space M of a dynamical system with N degrees of frecidom has a differentiable manifold structure, 
and the Lagrangian coordinates (gi, . . . , gjv) can be regarded as local coordinates on M. The latter becomes a 
Riemannian manifold once a proper metric is defined. For the sake of simplicity, let us consider systems of the form 
(1), so that the kinetic energy matrix is diagonal and the masses are all equal to one, i.e., aij = Sij. If we write 

gij=2[E-V{q)]dij ; (9) 

then Eq. (8) becomes 



0: 



j2Tdt = S j {gijq'-q^f'^ dt = 5 j ds , (10) 



so that the natural motions are the geodesies of M provided ds is the arc-length element, i.e., the metric on M is 
given by the tensor whose components are just the gij defined in Eq. (9). This metric is referred to as the Jacobi 
metric, and its arc-length element is 

ds^ = gijdqW = 2[E - V{q)] dt^ = 4[E - V{q)]^ dt" . (11) 

The geodesic equations written in the local coordinate frame {q^,. . . , q^) are (see Appendix A, Eq. (A34)) 

Di _ , pi dqi dq'' _ 

where D/ds is the covariant derivative along the curve 7(5) (see Appendix A, Eqs. (A27) and (A32)), 7 = dq/ds is the 
velocity vector of the geodesic and the T are the Christoffel symbols. Using the definition of the Christoffel symbols 
(see Appendix A, Eq. (A31)) and Eq. (9) it is straightforward to show that the Eqs. (12) become 



''Such a parameter is the time t here, but could also be the arc-length s in the following. 

^The accessible configuration space Me can then be seen as the union of all the "sub-configuration spaces" {g € M : V{q) = 
E — T} that one gets for all the possible values oiT, <T < E. 
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+ 



2{E-V) 



d{E - V) dq^ dq^ „■ d{E - V) dq'' dq"" 



ds ds 



dqj 



-Qkn 



ds ds 



= 



whence, using Eq. (11), Newton's equations are recovered, 

(fq^ 



df^ 



dV 
dqi ' 



(13) 



(14) 



Note that the Jaeobi metric is obtained by a conformal changx^ of the kinetic energy metric ciij sec Eq. (9) and 
Appendix A, §A3. In fact the general result for the Riemannian geometrization of natural Hamiltonian dynamics is 
the following 

Theorem. Given a dynamical system on a Riemannian manifold {M,a), i.e., a dynamical system whose Lagrangian 

is 

r. = laijq'q^ - V{q) , 
then it is always possible to find a conformal transformation of the metric, 



9^] 



such that the geodesies of (M, g) are the trajectories of the original dynamical system; this transformation is defined 
by 



<f{q) = log[E - Viq)] . 



The proof proceeds as above, using Eqs. (9)-(14) and simply replacing all the 6ij matrices with the kinetic energy 
matrix a^-; for details, see e.g. Ref. [11]. 



2. The Eisenhart metric 



We could try to consider the configuration spacetime M x R as an alternative ambient space for the geometrization 
of dynamics, with local coordinates (g" = t,q^,. . . ,q^), and to define a metric starting from Hamilton's principle 
S J £dt = 0. We could try to define a metric tensor by multiplying the Lagrangian (4) by 2{dq^Y: 

ds^ = 2C {dq°f = {gc)^, dq^dq"" = aij dq'dq^ - 2V{q) {dq°f , (15) 

where and v run from to and i and j run from 1 to N; however, one can easily verify that the geodesies of the 
manifold (M x R, gc) are then not the natural motions of the systems, since the Lagrangian is not an integral of the 
motion. 

However, we can consider an ambient space with an extra dimension, M x R^, with local coordinates 
((f ,q^ , . . . ,q^ , . . . ^q'^ ,q^^^). This space can be endowed with a non-degenerate pseudo-Riemannian metric (see 
Appendix A, §A lb)), first introduced by Eisenhart [52], whose arc- length is 

c^s" = i3E)^, dq'^dq'^ = a,, dq'dq^ - 2V{q){dq'>f + 2 dq''dq^+^ , (16) 

where /j, and v run from to + 1 and i and j run from 1 to N, and which, from now on, will be referred to as the 
Eisenhart metric, and whose metric tensor will be denoted as gs- The relation between the geodesies of this manifold 
and the natural motions of the dynamical system is contained in the following [53] 

Theorem (Eisenheirt) . The natural motions of a Hamiltonian dynamical system are obtained as the canonical 
projection of the geodesies of (M x R^, gs) on the configuration space-time, n : M x 'R? i— > M x R. Among the totality 
of geodesies, only those whose arclengths are positive- definite and are given by 

ds"^ = cldt^ (17) 

correspond to natural motions; the condition (17) can be equivalently cast in the following integral form as a condition 
on the extra- coordinate q^~^^: 
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-/' 



CdT 



(18) 



where ci and C2 are given real constants. Conversely, given a point P G M x R belonging to a trajectory of the system, 
and given two constants ci and C2, the point P' = ■jt~^{P) e M x R^, with q^~^^ given by (18), describes a geodesic 

curve in (M x R^, g£;) such thai ds"^ = c^^df^ . 

For the full proof, see Ref. [53]. Since the constant Ci is arbitrary, we wiU always set = 1 in order that ds"^ = dt^ 
on the physical geodesies. 

From Eq. (16) follows that the explicit table of the components of the Eisenhart metric is 



/ -2V^(?) 
an 



9e 



\ 



oni 




l\ 
aiN 



onn 
/ 



(19) 



where is the kinetic energy metric. The non- vanishing Christoffel symbols, in the case = 5ij , are only 

rj,o = -r^' = diV , (20) 



so that the geodesic equations (12) read 



'2„0 



, 



2„i 



ds"^ 



ds 



(21a) 
(21b) 



using ds = dt one obtains 



'2_0 



pAT+l dq" dq' _ 
°' ds ds - ^ 

= , 

dV 
dqi 



df^ 



df^ 



(21c) 

(22a) 

(22b) 



dt^ 



dC 
"dt 



(22c) 



Eq. (22a) only states that q^ = t, The A'' equations (22b) are Newton's equations, and Eq. (22c) is the differential 

version of Eq. (18). 

The fact that in the framework of the Eisenhart metric the dynamics can be geometrized with an afhne parametriza- 
tion of the arc- length, i.e., ds = dt, will be extremely useful in the following, together with the remarkably simple 
curvature properties of the Eisenhart metric (see §111 C). 



B. Curvature and stability 

The gcometrization of the dynamics is a natural framework for the study of the stability of the trajectories of a 
dynamical system, for it links the latter with the stability of the geodesies; the latter is completely determined by the 
curvature of the manifold, as shown below. 

Studying the stability of the dynamics means determining the evohition of perturbations of a given trajectory. This 
implies that one should follow the evolution of the linearized (tangent) flow along the reference trajectory. For a 
Newtonian system, writing the perturbed trajectory as 
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(23) 



substituting this expression in the equations of motion 

and retaining terms up to first order in the ^'s, one finds that the perturbation obeys the so-called tangent dynamics 
equation which reads as 

This equation should be solved together with the dynamics in order to determine the stability or instability of the 
trajectory: when the norm of the perturbations grows exponentially, the trajectory is unstable, otherwise it is stable. 

Let us now translate the stability problem into geometric language. By writing, in close analogy to what has been 
done above in the case of dynamical systems, a perturbed geodesic as 

q\s) = q\s) + J\s) , (26) 

and then inserting this expression in the equation for the geodesies (12), one finds that the evolution of the perturbation 
vector J is given by the following equation: 

DV* i dq^ ^dq^ 

-^+R.^,—J—-0, (27) 

where R^j^i are the components of the Riemann curvature tensor (see Appendix A, Eq. (A39)). Equation (27) is 
referred to as the Jacobi equation, and the tangent vector field J as the Jacobi field. This equation was first studied 

by Lcvi-Civita and is also often referred to as the equation of Jacobi and Lcvi-Civita. For a derivation we refer to 
Appendix A, §A4, where it is also shown that one can always assume that J is orthogonal to the velocity vector along 
the geodesic, 7, i.e., 

(J,7) = 0, (28) 

where (•, •) stands for the scalar product induced by the metric (see Appendix A, Eq. (A20)). The remarkable fact 
is that the evolution of J — and then the stability or instability of the geodesic — is completely determined by 
the curvature of the manifold. Therefore, if the metric is induced by a physical system, as in the case of Jacobi or 
Eisenhart metrics, such an equation links the stability or instability of the trajectories to the curvature of the ambient 
manifold. 

The subject of the next sections is precisely to exploit such a link in order to describe and understand the behaviour 
of those physical systems whose trajectories are mainly unstable. 

However, before that, we have to give explicit expressions for the curvature of the mechanical manifolds, i.e., of 
those manifolds whose Riemannian structure is induced by the dynamics via the Jacobi or the Eisenhart metric. 



C. Curvature of the mechanical manifolds 



We already observed that the Jacobi metric is a conformal deformation of the kinetic-energy metric, whose compo- 
nents are given by the kinetic energy matrix aij . In the case of systems whose kinetic energy matrix is diagonal, this 
means that the Jacobi metric is conformally flat (see Appendix A, §A3). This greatly simplifies the computation of 
curvatures. It is convenient to define then a symmetric tensor C whose components are [11] 



N-2 

4(E-y)2 



2{E - V)didjV + SdiVdjV 



(29) 



where V is the potential, E is the energy, and V and | • | stand for the Euclidean gradient and norm, respectively. 
The curvature of {ME,gj) can be expressed through C. In fact, the components of the Riemann tensor are 



1 



■ijkm 



N-2 



(30) 
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By contraction of the first and third indices, we obtain the Ricci tensor, whose components are (see Appendix A, Eq. 
(A46)) 

= AiE-Vr f^^^ ~ ^^^'^'^ ^ mVdjV] + j^^^ [^{E - V)AV - (TV - 4)| VFfl , (31) 
and by a further contraction we obtain the scalar curvature (see Appendix A, Eq. (A48)) 

The curvature properties of the Eisenhart metric are much simpler than those of the Jacobi metric, and this is 
obviously a great advantage from a computational point of view. The only non- vanishing components of the curvature 
tensor are 



RoiOj = didjV , (33) 
hence the Ricci tensor has only one nonzero component 

Roo = AV , (34) 

and the scalar curvature is identically vanishing, 

n = . (35) 

To summarize, we have shown that the dynamical trajectories of a Hamiltonian system of the form (1) can be seen 
as geodesies of the configuration space, or of an enlargement of it, once a suitable metric is defined^. The general 
relationship which holds between dynamical and geometrical quantities regardless of the precise choice of the metric 
can be sketched as follows: 

dynamics geometry 



(time) t ~ s (arc-length) , . 

(potential energy) V ^ g (metric) ^ ' 

(forces) dV ^ T (Christoffel symbols) 

("curvature" of the potential) d^V,{dV)'^ ~ R (curvature of the manifold) 

In the case of the Eisenhart metric, all these relations are extremely simple (maybe as simple as possible). In fact the 
physical time t can be chosen as equal to the arc-length s, the metric tensor cje contains only the potential energy V, 
the non- vanishing Christoffel symbols F are equal to the forces dV, and the components of the Riemann curvature 
tensor R contain only the second derivatives of the potential energy, didjV. 

We have also shown that the stability of the dynamical trajectories can be mapped onto the stability of the geodesies, 
which is completely determined by the curvature of the manifold. We will show in §IVA3 that, in the case of the 
Eisenhart metric, as a consequence of its remarkably simple properties, also the relationship between the stability of 
the trajectories and the stability of the geodesies becomes as simple as possible, i.e., the Jacobi equation (27) becomes 
identical to tangent dynamics equation (25). 



IV. GEOMETRY AND CHAOS 



The purpose of the present section is to describe in some detail how it is possibile, using the Jacobi equation as the 
main tool, to reach a twofold objective: first, to obtain a deeper understanding of the origin of chaos in Hamiltonian 

systems, and second, to obtain quantitative informations on the "strength" of chaos in these systems. Some basic 
concepts about Hamiltonian chaos and the definition of Lyapunov exponents are summarized in Appendix C. 



As already stated at the beginning of this Section, there are many other possible choices for the ambient manifold and it 
metric: some other possible choices are described in Ref. [33]. 



11 



A. Geometric approach to chaotic dynamics 



A physical theory should provide a conceptual framework for modeling and understanding — at least at a qualitative 
level the observed features of the system which is the object of the theory, and should also have a predictive content, 
i.e. should provide quantitative tools apt to compute, at least approximately, the outcomes of the experiments (no 
matter if it concerns laboratory experiments or numerical experiments performed on a computer). According to these 
requirements, a satisfactory theory of deterministic chaos is certainly still lacking. In fact in both aspects the current 
theoretical approaches to chaos have some problems, especially if we consider the case of conservative flows, i.e., of 
the dynamics of conservative systems of ordinary differential equations. 

To explain the origin of chaos in conservative dynamics one usually invokes the existence of invariant hyperbolic 
sets — or horseshoes — in phase space, like those generated by homoclinic intersections of perturbed separatrices (see 
Appendix C). In order to quantify the degree of instability of a trajectory or of a system we must instead resort to 
the notion of Lyapunov exponents. The Lyapunov exponents are asymptotic quantities and their relation with local 
properties of phase space, like horseshoes, is far from evident; nonetheless they provide the natural measure of the 
degree of chaos, measuring the typical time scales over which a trajectory looses the memory of its initial conditions. 
A rigorous definition of the existence of chaotic regions in the phase space of a system, based on the detection of 
horseshoes, does not provide any quantitative tool to measure chaos; on the other hand, Lyapunov exponents allow 
a very precise measure of chaos but give no information at all on the origin of such a chaotic behaviour. From 
a conceptual point of view this situation is far from being satisfactory, not to speak of the fact that the practical 
application of the methods to search horseshoes becomes extremely difficult as the number of degrees of freedom 
is large [54] . From the predictive point of view the situation is even worse, for no analytic method at all exists to 
compute Lyapunov exponents, at least in the case of flows of physical relevance. It is worth noticing that in a recent 
paper [55], Gozzi and Renter have shown that one could build, in principle, a field-theoretic framework to compute 
Lyapunov exponents, but the practical application of such methods is still unclear. Needless to say, all the tools 
belonging to canonical perturbation theory, which have undergone remarkable developments in the last years [56], 
can hardly be used to compute quantities like Lyapunov exponents since in this framework one can only describe the 
regular, i.e., nonchaotic, features of phase space. 

The geometric approach to dynamical instability allows a unification of the method to measure chaos with the 
explanation of its origin. In fact the evolution of the field J given by the Jacobi equation (27) contains all the 
information needed to compute Lyapunov exponents, and makes us also recognize in the curvature properties of the 
ambient manifold the origin of chaotic dynamics. 

Obviously, also this approach is far from being free of problems. For instance, the only case in which it is possible 
to rigorously prove that some definite curvature properties imply chaos in the geodesic flow, is the case of compact 
manifolds whose curvature is everywhere negative. In this case every point of the manifold is hyperbolic: in a sense 
this is the opposite limit to the integrable case. Though abstract and unphysical, such systems can help intuition. In 
a geodesic flow on a compact negatively-curved manifold, the negative curvature forces nearby geodesies to separate 
exponentially, while the compactness ensures that such a separation does not reduce to a trivial "explosion" of the 
system and obliges the geodesies to fold. The joint action of stretching and folding is the essential ingredient of chaos. 

Krylov tried to apply this framework to explain the origin of mixing in physical dynamical systems. Unfortunately 
for many systems in which chaos is detected the curvatures are found mainly positive, and there are examples, for 
instance the Hcnon-Heiles system — see Eq. (42) — geometrized with the Jacobi metric and the Fermi-Pasta-Ulam 
model — see Eq. (83a) — geometrized with the Eisenhart metric, where curvatures are always positive even in the 
presence of fully developed chaos. Hence even positive curvature must be able to produce chaos. 

Only recently an example has been found of a compact surface with positive curvature, where the presence of chaotic 
regions coexisting with regular ones can be rigorously proved [57], and this provides mathematical support for the 
available numerical evidence that negative curvature is not necessary at all to have chaos in a geodesic flow [22,24,31]. 
What then is the crucial feature of the curvature which is required to produce chaos? There is not yet a deflnite 
answer at least on rigorous grounds to this question. Nevertheless it is sure that, if positive, curvature must be 
non-constant in order to originate instability, and we shall see that the curvature fluctuations along the geodesies can 
be responsible for the insurgence of an instability through a mechanism very close to parametric instability. 

The advantages of the geometric approach to chaos are not only conceptual: also on predictive grounds this 
framework proves very useful. For, starting from the Jacobi equation, it is possible to obtain an effective stability 
equation which allows one to obtain an analytic estimate of the largest Lyapunov exponent in the thermodynamic 
limit [25,31]. Such an estimate turns out to be in very good agreement with the results of numerical simulations for 
a number of systems (see §IVB). In order to understand the derivation of such an efl[ective stability equation, let us 
investigate in greater detail the relation between stability and curvature which was introduced in the last section. 
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1. Geometric origin of Hamiltonian chaos 



Let us consider an A/^-dimensional Riemannian (or pseudo-Riemannian) manifold (M, g) and a local coordinate 
frame with coordinates (q^ , . . . , g^). 

We already observed that the evolution of the Jacobi field J, which contains the whole information on the stability of 
the geodesic flow, is completely determined by the curvature tensor R through the Jacobi equation (27). Unfortunately 
the number of independent components of the tensor R is 0{N'^) — even if this number can be considerably reduced 
by symmetry considerations — so that Eq. (27) becomes rather untractable already at fairly small dimensions. 

Nevertheless there is a particular case in which the Jacobi equation has a remarkably simple form: the case of 
isotropic manifolds (see Appendix A, §A3), where Eq. (27) becomes 

—^+KT = 0, (37) 

where K is the constant sectional curvature of the manifold (see Appendix A, Eq. (A42)). Choosing a geodesic frame, 
i.e. an orthonormal frame transported along the geodesic, covariant derivatives become ordinary derivatives, i.e., 
D/ds = d/ds, so that the solution of Eq. (37), with initial condition J(0) = and dJ{0)/ds = w{0), is 

' ^sm(^VKs) {K>0); 

J{s) = { sw{s) {K = 0) ; (38) 

^smh{^/^s) {K<0) . 

The geodesic flow is unstable only if K < 0, and in this case the instability exponent is just y/—K. 

As long as the curvatures are negative, the geodesic flow is unstable even if the manifold is no longer isotropic, 
and by means of the so-called "comparison theorems" (mainly Ranch's theorem, see e.g. [58]) it is possible to prove 
that the instability exponent is greater or equal to (— maxM(-ft^))^^^ [7]. On the contrary, no exact results of general 
validity have yet been found for the dynamics of geodesic flows on manifolds whose curvature is neither constant nor 
everywhere negative. 

Equation (37) is valid only if K is constant. Nevertheless in the case in which dimM = 2 (surfaces), the Jacobi 
equation again written in a geodesic reference frame for the sake of simplicity — takes a form very close to that 
for isotropic manifolds, 

fl+K{s)J = 0, (39) 

where 

K{s) = ^n{s) (40) 

and. contrary to Eq. (37), it is no longer a constant. With TZ{s) we denote the scalar curvature of the manifold at the 
point P = j{s) (see Appendix A, Eq. (A45)). If the geodesies are unstable. Equation (39) has exponentially growing 
solutions. As far as we know, the solutions of (39) can exhibit an exponentially growing envelope in two cases: 

a. the curvature K{s) takes negative values; 

b. the curvature K{s), though mainly or even exclusively positive, fluctuates in such a way that it triggers a sort 
of parametric instability mechanism. 

In the first case, the mechanism that is at the origin of the instability of the geodesies is the one usually considered 
in ergodic theory [7]. But in the second case a new mechanism of instability, which does not require the presence of 
negatively curved regions on the manifold, shows up: the fluctuations of the curvature along the geodesic make the 
geodesic unstable. 

Let us now turn to physics, i.e., to the case of a mechanical manifold: in the case of the Jacobi metric with N = 2 
the scalar curvature written in standard (Lagrangian) coordinates reads as 

{E - y)3 ^ {E- vy ' ^ ' 
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where V and A stand respectively for the euchdean gradient and Laplacian operators. Hence we can have 7?. < only 
if AV < 0, i.e., for stable physical potentials, when the potential has inflection points. In these cases Krylov's idea 
can work — even if in the high-dimensional case this becomes very complicated — and we may have dynamical chaos 
induced by negative curvatures of the manifold. Indeed Krylov was mainly concerned with weakly non-ideal gases, or 
in general dilute systems, where for typical interatomic interactions AV < so that the curvatures can be negative 
(see Krylov's PhD thesis in Ref. [5]). 

We will now show one example in which, though chaos is present, curvatures are positive. Let us consider the 
Henon-Heiles model [59], whose Hamiltonian is 

^ = 1{pI+pI) + 1 + y') + ^'y - ly' ■ (42) 

This model was introduced in an astrophysical framework to study the motion of a star in an axially symmetric 

galaxy, but it can also bo regarded as a model of a triatomic molecule (after having used translational symmetry to 
eliminate the center-of-mass coordinate) [60] . The Henon-Hciles model is a cornerstone in the study of Hamiltonian 
chaos: it was the first physical model for which chaos was found and where a transition from a mainly regular to 
a mainly chaotic phase space was identified under a variation of the energy. In this model, Eq. (39) is exact, but 
TZ > everywhere. Hence chaos in this system cannot come from any negative curvature in the associated mechanical 
(Jacobi) manifold. As we shall sec later on (see e.g. §IVB), the absence of negative curvatures in the associated 
mechanical manifolds is not a peculiarity of this model, for it is shared with many systems of interest for field theory 
and condensed matter physics which have chaotic trajectories. In particular, all the systems that in the low-energy 
limit behave as a collection of harmonic oscillators do belong to this class. 

In these cases the second of the previously discussed instability mechanisms, the one mentioned in item b, comes 
in: curvature fluctuations may induce chaos through parametric instability. The latter is a well-known feature of 
differential equations whose parameters arc time-dependent. The classical example (see e.g. Arnol'd's book [2]) is the 
mathematical swing, i.e. a pendulum, initially at rest, whose length is modulated in time. If the modulation contains 
frequencies resonating with the free pendulum's fundamental frequency, the st able equilibrimn position gets unstable 
and the swing starts to oscillate with growing amplitude. In Eq. (39), \/K{s) and s play the role of a frequency and 
of time, respectively, so that this equation can be thought of as the equation of motion of a harmonic oscillator with 
time-dependent frequency, often referred to as a (generalized) Hill's equation [61]. By expanding K{s) in a Fourier 
series we get 

oo 

K{s) = Ko + [a„ cos{nujs) + bn sin(nws)] , (43) 

n=l 

where ui = 2'it/L and L is the length of the geodesic. The presence of resonances between the average frequency a/Kq 
and the frequency in some term in the expansion (43) eventually forces an exponential growth of the solutions of the 
equation. In the simplest case, in which only one coefficient of the series (43), say ai, is non- vanishing, the equation is 
called the Mathieu equation and it is possible to compute analytically both the bounds of the instability regions in the 
parameter space and the actual value of the instability exponent [61]. At variance with the Mathieu case, in the general 
case, where a large number of coefficients of the Fourier decomposition of K{s) is nonzero, it is much more difficult to 
do something similar. Hence there is not yet any rigorous proof of the fact that this kind of parametric instability is 
the mechanism that produces chaos in Hamiltonian dynamical systems — in the two-dcgrces-of-freedom case or in the 
general case — and this still remains a conjecture. Nevertheless such a conjecture is strongly supported by at least 
two facts. First of all, in recent papers [30,27] it has been shown that the solutions of the Jacobi equation (39) for the 
Henon-Heiles model and for a model of quartic coupled oscillators show an oscillatory behaviour with an exponentially 
growing envelope — which is precisely what one expects from parametric instability — in the chaotic regions, while 
the oscillations are bounded in the regular regions. Second, also in high- dimensional flows the components of the 
Jacobi field J oscillate with an exponentially growing amplitude as long as the system is non-integrable, whereas they 
exhibit only boimded oscillations for integrable systems. Moreover, in the high-dimensional case (i.e., for systems 
with a large number of degrees of freedom) it is possible to establish a quantitative link between the largest Lyapunov 
exponent and the curvature fluctuations. In fact, as we shall see in the following, in the high-dimensional case it 
is possibile to write down, under suitable approximations, an effective stability equation which looks very similar to 
Eq. (39), but where the squared frequency K(s) is a stochastic process, and, through this equation, it is possible to 
give an analytical estimate of the largest Lyapunov exponent. Since, from now on, we are going to consider only the 
largest Lyapunov exponent, the latter will be referred to as just the Lyapunov exponent. 
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2. Effective stability equation in the high- dimensional case 



Let us now study the problem of the stabihty of the geodesies in manifolds whose dimension N is large: according 
to the correspondence between geometry and dynamics introduced in Sec. Ill, we are considering a system with a 
large number A'^ of degrees of freedom. 

Our starting point is the Jacobi equation (27). Our aim is to derive from it an effective stability equation which no 
longer depends on the dynamics, i.e., on the evolution of the particular geodesic that we are following, but only on 
the average curvature properties of the manifold. To do that, we need some assumptions and approximations which 
are not valid in general but which are very reasonable in the case of large- A'' mechanical manifolds. For the sake of 
clarity, we first summarize the assumptions and approximations leading to our final result, and later on we discuss 
them more thoroughly. Further details can be found in the papers where this approach was originally put forward 
[25,31]. 

0. We assume that the evolution of a generic geodesic is chaotic. This assumption is reasonable in the case of a 

manifold whose geodesies are the trajectories of a generic Hamiltonian system with a large number of degrees of 
freedom iV, for in this case the overwhelming majority of the trajectories will be chaotic. This bears a certain 
similarity to Gallavotti and Cohen's "Chaotic Hypothesis" [62]. 

1. We assume that the manifold is quasi-isotropic. Loosely speaking, this assumption means that the manifold can 

be regarded somehow as a locally deformed constant-curvature manifold. However, we will give this assumption 
a precise formulation later, in Eqs. (49). This approximation allows us to get rid of the dependence of the Jacobi 
equation (27) on the full Riemann curvature tensor by replacing it with an effective sectional curvature /C(s) 
along the geodesic; moreover, the Jacobi equation becomes diagonal. 

2. To get rid of the dependence of the effective sectional curvature /C(s) on the dynamics, i.e., on the evolution of 
the geodesic, we model /C(s) with a stochastic process. This assumption is motivated by Assumption above. 
Moreover, as long as we consider a high-dimensional mechanical manifold associated to a Hamiltonian flow 
with A'^ degrees of freedom and we are eventually interested in taking the thermodynamic limit N oo, the 
sectional curvature is formed by adding up many independent terms, so that invoking a central-limit-theorem-like 
argument, /C(s) is expected to behave, in first approximation, as a gaussian stochastic process. 

3. We assume that the statistics of the effective sectional curvature JC is the same as that of the Ricci curvature 
Kji, which is a suitably averaged sectional curvature (see Appendix A, Eq. (A44)). Such an assumption is 
consistent with Assumption 1 above, for in a constant-curvature manifold the sectional curvature equals the 

Ricci curvature times a constant, and allows us to compute the mean and the variance of the stochastic process 
introduced in Assumption 2 in terms of the average and the variance of Kn along a generic geodesic. 

4. The last step, which completely decouples the problem of the stability of the geodesies from the evolution of the 
geodesies themselves, consists in replacing the (proper) time averages of the Ricci curvature with static averages 
computed with a suitable probability measure /i. If the manifold is a mechanical manifold, the natural choice 
for II is the microcanonical measure. Again this assumption is reasonable if Assumption is valid. 

After these steps, we end up with an effective stability equation which no longer depends on the evolution of the 
geodesies, but only on the average and fluctuations of the Ricci curvature of the manifold. 

Let us now discuss more thoroughly the above-sketched procedure. For that, it is convenient to introduce the Weyl 
projective tensor W, whose components are given by [63] 

W%i = R'^u - {RjA - RjkSi) , (44) 

where Rij = R^l^j are the components of the Ricci curvature tensor (see Appendix A, Eq. A46). Weyl's projective 
tensor measures the deviation from isotropy of a given manifold, since it vanishes identically for an isotropic manifold. 
Then we can reformulate the Jacobi equation (27) in the following form [31]: 

D^J' . 1 dq^dq" 1 dq^ ,dq^ dq^ ,dq^_ 

-d^ + N^^^'^^-^ -N^^^"^-^ (^^^ 

For an isotropic manifold the third term in Eq. (45) vanishes because Rjk = K gjk (see Appendix A, Eq. A50) so that 
RjkQ^ J'^ = K (7, J), and (7, J) = (see Eq. 28). Thus, for an isotropic manifold Eq. (45) collapses to Eq. (37), in fact 
the second term is nothing but K J^. When the manifold is not isotropic, we see that Eq. (45) retains the structure 
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of Eq. (37) up to its second term, since the coefficient of J* is still a scalar. This coefficient has now the meaning of a 
sectional ciirvature averaged, at any given point, over the — 1 independent directions orthogonal to 7, the velocity 
vector of the geodesic. However, such a mean sectional curvature is no longer constant along the geodesic 7(5), and 
is nothing but that the Ricci curvature Kr divided hy N — 1 (see Appendix A, Eq. A47). The fourth term of (45) 
accounts for the local degree of anisotropy of the ambient manifold. 

Let us now consider a geodesic frame: in this case Eq. (45) can be rewritten as 

^ + fcfl(s) r - r] (s) + w] (s) = (46) 
where, to ease the notation, we have put 

I. ( \ 1 D dq' dq'' 



Being a scalar quantity, the value of Ur is independent of the coordinate system. Now let us formulate our Assumption 
1, namely, that the manifold is quasi-isotropic, in a more precise way. To do that, we recall that (see Appendix A, 
Eqs. (A49) and (A50)) if and only if the manifold is isotropic, i.e., has constant curvature, the Riemann curvature 
tensor and the Ricci tensor can be written in the remarkably simple forms, i.e., 

Rijki = K {giugji - gugjk) , (48a) 

and 

Rij=Kgij, (48b) 

where if is a scalar constant, the sectional curvature of the manifold. The precise formulation of Assumption 1 is 
now that along a generic geodesic the Riemann curvature tensor and the Ricci tensor retain the same functional form 
as in the case (48), i.e., that 

Rijki ~ ^(s) {gikgji - gugjk) , (49a) 

and 

Rij « /C(s) g^J , (49b) 

where /C(s), which is no longer a constant, is an effective sectional curvature. In the general case we are not able to 
give a rigorous explicit expression for /C(.s), because the functional dependence postulated in Eqs. (49) holds true only 
for constant-curvature manifolds. However, the effective curvature /C(s) is expected to be essentially the sectional 
curvature K{j, J) (see Appendix A, Eq. (A44)) measured along the geodesic in the directions of the velocity vector 
j — dq/ds and of the Jacobi vector J. 

Combining Eqs. (47b) and (49b), and recalling that the vector J is orthogonal to the velocity of the geodesic, i.e., 

gij^J-' = 0, we find that the third term in Eq. (46), —rjJ\ vanishes as in the isotropic case. Now we combine Eqs. 
(44) and (49a) to obtain 

W),i ^ JC{s){5]gui - 5\gkj) - j^iRji^l - RjkS}) (50) 
so that Eq. (47c) can be rewritten as 

w] ^ lC{s)5] - kn{s)5) - /C(s)^^5.,-^ , (51) 



16 



where we have used the definition of kR given in Eq. (47a) and the approximation (49b) for the Ricci tensor. Let us 
now insert Eq. (51) into Eq. (46): the last term of Eq. (51) vanishes after having been multipUed by and summed 
over j, beeause J and dq/ds are orthogonal, and the term kii{s)J^ is canceled the term — fci{(s)J' coming from Eq. 
(51), so that Eq. (46) is finally rewritten as 



d^P 
ds"^ 



/C(s) r = 0. (52) 



Equation (52) is now diagonal. However, in order to use it, we should know the values of JC{s) along the geodesic. 
Here, Assumptions 2 and 3 come into play: we replace /C(s) with a stochastic gaussian process, and we assume that 
its probability distribution is the same as that of the Ricci curvature, 

P(/C) ~ ViKR) . (53) 

Such an assumption is consistent with our Assumption 1, because for an isotropic manifold the sectional curvature is 
identically equal to the Ricci curvature divided by — 1, so that, if the manifold is quasi-isotropic, it is natural to 
assume that the probability distributions of the sectional curvature and of the Ricci curvature are similar. Moreover, 
such an assumption is also the only easy one, because we are able to compute, under some further assumptions, the 
probability distribution of Kr, but we do not know anything about /C. 

To be consistent with the definition of the sectional and the Ricci curvatures (see Appendix A, Eq. (A51)), the 
following relations are assumed to hold for the first two cumulants of (53): 

■j^{Kr{s))s ^ (kRis)), , (54a) 



ms) - tCr), ^{[Kr{s) - {Kr),]'), ^ {6'kR)s , (54b) 

where (•)<, stands for a proper-time average along a geodesic 7(s). In general, the probability distributions (53) will 
not be Gaussian, i.e., other cumulants in addition to the first two will be nonvanishing. However, we already observed 
that since for a large system Kr is obtained by summing a large number of independent components, it is reasonable 
to assume that a sort of central limit theorem holds and that a gaussian approximation is sufficient. 
Our approximation for the effective sectional curvature /C(s) is then the stochastic process 

/C(s) ~ {kR{s))s + (5'fci?)P vis) , (55) 

where 77(5) is a random gaussian process with zero mean and unit variance. 

Finally, in order to completely decouple the stability equation from the dynamics, we use Assumption 4 and we 
replace time averages with static averages computed with a suitable measure If the manifold is a mechanical 
manifold the geodesies are the natural motions of the systems, and a natural choice for is then the microcanonical 
ensemble, so that Eq. (55) becomes 

]C{s)^{kR{s))^ + {S'kR)l/'rjis) . (56) 
Our final effective stability equation is then 

+ (fcfl)^ ^ + {5^kR)]!^ vis) V = , (57) 

where ■ip stands for any of the components of J, since all of them now obey the same effective equation of motion. 

Equation (57) implies that, if the manifold is a mechanical manifold, the growth-rate of tp gives the dynamical 
instability exponent in our Riemannian framework. Equation (57) is a scalar equation which, independently of the 
knowledge of the dynamics, provides a measure of the degree of instability of the dynamics itself through the behavior 
of ipis). The peculiar properties of a given Hamiltonian system enter Eq. (57) only through the global geometric 
properties {kR)^ and {S^kR}^ of the ambient Riemannian manifold (whose geodesies are natural motions) and are 
sufficient, as long as our Assumptions 1-4 hold, to determine the average degree of chaoticity of the dynamics. 
Moreover, (fc/j)^ and {5'^kR) ^ are microcanonical averages, so that they arc functions of the energy E of the system, 
or of the energy per degree of freedom e = E/N which is the relevant parameter as iV ^ 00. Thus from (57) we can 
obtain the energy dependence of the geometric instability exponent. 
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Within the vahdity of our Assumptions 1-4, transforming the Jacobi equation (27) into Eq. (57), its original 
complexity of the Jacobi equation has been considerably reduced: from a tensor equation we have obtained an 
effective scalar equation formally representing the equation of motion of a stochastic oscillator. In fact, Eq. (57), with 
a self-evident notation, is of the form 

^ + k{s)i; = (58) 

where k{s), the squared frequency, is a gaussian stochastic process. 

Moreover, such an equation admits a very suggestive geometric interpretation, since it is a scalar equation, i.e., it 
is formally the .Jacobi equation on a 2-dimensional manifold whose Gaussian curvature is given, along a geodesic, by 
the random process k{s), which can be regarded as an "effective" low-dimensional manifold approximatng the "true" 
high-dimensional manifold where the dynamics of the geodesic flow takes place. This is the real geometrical content of 
our quasi-isotropy hypothesis. Hence the average global curvature properties (fc/j)^ and ((5^fc/f)^, in addition to being 
the ingredients for a geometric computation of the instability exponent, convey also information on the geometric 
structure of this eflfective manifold. Thus we expect that it will be possible to gain some insight in the global properties 
of the dynamics by simply studying the behaviour of these average curvature properties as the energy is varied. 



3. A geometric formula for the Lyapunov exponent 

Let us now study the properties of the solutions of Eq. (58) in order to obtain an analytic estimate for the Lyapunov 
exponent. The derivation of the stochastic oscillator equation does not depend on a particular choice of the metric; 
within the approximations discussed above. Equation (58) holds regardless of the choice of the metric. However, to 
make explicit the connection between the solutions of Eq. (58) and the stability of a dynamical system, one has to 
choose a particular metric; in the case of Hamiltonian systems of the form (1), the choice of the Eisenhart metric is 
the simplest one. 

For this reason, we shall from now on restrict ourselves to standard Hamiltonian systems with a diagonal kinetic 
energy matrix, i.e., Uij = Sij, choosing as ambient manifold for the geometrization of the dynamics the enlarged 
configuration space-time equipped with the Eisenhart metric (16). The case of the Jacobi metric is discussed in Ref. 
[31]. 

We will proceed as follows, (i) We will show that in the present case the Jacobi equation (27) is equal to the 
tangent dynamics equation (25). (ii) We will replace the arc-length s with the time t and we will explicitly compute 
the average and the fluctuations of the Ricci curvature along a geodesies in terms of dynamical observables, so that 
the (static) probability distribution of the stochastic process k{t) which models the effective sectional curvature is 
defined, {in) We will give an estimate for the time correlation function of the process k{t). {iv) We will solve the 
stochastic oscillator equation, obtaining an analytical formula for the Lyapunov exponent. 

Let us now consider item (i). As a consequence of the simple structure of the curvature tensor for the Eisenhart 
metric (sec §111 C), the Jacobi equation (27) takes the form (we recall that the manifold has now dimension A'' -|- 2; 
all the indices k, . . . run from 1 to A^) 

D^JO dq^ dq^ dq' dq^ _ 



D'r [dq'\\,.j.. dq" ,dq^ , dq^ ,dq'_ 

dq° 

-J^^=0, (59c) 



, j.N+ldq' rodq' , r,N+ldq' dq^ 
ds^ + ^o^' ds ^^^-'O rfs ' 



where, for the sake of clarity, we have written out Eq. (59) separately for the 0, the i = 1,. . . ,N, and the A^ -|- 1 
components, respectively. As F?- = (see Eq. 20) we obtain, from the definition of covariant derivative (see Appendix 
A, Eq. A32), DJ^/ds = dJ^/ds, and, as R^i^,^ = (see §niC), we find that Eq. (59a) becomes 

^=0, (60) 
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so that j" does not accelerate and, without loss of generality, we can set ^4- = J'^(O) = 0. Combining the latter 
result with the definition of covariant derivative we obtain 

DJ^ dJ' , dq^ 



^Ik^J" (61) 



and using dq^ /ds = we finally get 



ds ds ds 



(62) 



ds'^ ds"^ 

so that Eq. (59b) gives, for the projection in configuration space of the separation vector, 

ds"^ dqidqk \ ds J 

Equation (59c) describes the evolution of J^+^, which, however, does not contribute to the norm of J because 
9N+1N+1 = 0, so we can disregard it. 

Along the physical geodesies of gs, ds^ = {dq^Y — '^^^ ■> so that Eq. (63) is exactly the usual tangent dynamics 
equation (25), provided the identification ^ = J is made. This clarifies then the relationship between the geometric 
description of the instability of a geodesic flow and the conventional description of dynamical instability. We stress 
that from a formal viewpoint this is a peculiarity of the Eisenhart metric; nevertheless the physical content of this 
result is valid independently of the metric used, as long as the identification between trajectories and geodesies holds 
true. For, in recent papers [30,27] it has been found that using the .Jacobi metric the solutions of the .Jacobi equation 
and of those of the tangent dynamics equation — which in this case are two well-distinct equations — look strikingly 
similar. 

We now turn to item {ii). The Ricci curvature is obtained saturating the Ricci tensor with the components of the 
velocity vector dq/ds (see Appendix A, Eq. (A47)). In the present case, the only non- vanishing component of the 
Ricci tensor is i?oo = (see Eq. (34)), so that the dynamical observable which corresponds to the Ricci curvature 
along a geodesic depends only on the coordinates and not on the velocities and reads 

KR{q) = AV , (64) 

where we have used that, along a physical geodesic. {dq^Y = dt^ = ds^ ■ Using again this result we replace the 
arc-length s along the geodesic with the physical time t, and the stochastic oscillator equation (58) can be written 



dt^ 

where mean and variance of k{t) are given by 



$ + fc(i)V = 0, (65) 



fco = {kR)^ = ^(Ay)^ , (66a) 



al ^ (S'kn), = 1 {{{AVf), - {AV)l) . (66b) 

Since we are considering systems with large N — eventually taking the limit N ^ oo — we replaced TV — 1 with A'^ 
inEqs. (66). 

We consider now item (iii). The process k{t) is not completely defined unless its time correlation function, 

Tk{tut2) = {k{ti)k{t2)) - {k{h)){k{t2)) , (67) 

is given. The simplest choice is to assume that k{t) is a stationary and 5-correlated process, so that 

rk{ti,t2)=Tk{\t2-h\)=Tk{t)=Tal5{t) , (68) 

where r is the characteristic correlation time scale of the process. 

Before we can actually solve Eq. (65), we have then to give an explicit expression for r. To do that, first we will 
show how two independent characteristic correlation time scales, which will be referred to as ti and T2, respectively, 
can be defined, then we will estimate r by combining these two time scales. 
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A first time scale, which we will refer to as ri, is associated to the time needed to cover the average distance between 
two successive conjugate points along a geodesic. Conjugate points [58] are the points where the Jacobi vector field 
vanishes. As long as the curvature is positive and its fluctuations are small, compared to the average, two nearby 
geodesies will remain close to each other until a conjugate point is reached. At each crossing of a conjugate point 
the Jacobi vector field increases as if the geodesies were kicked (this is what happens when parametric instability is 
active). Thus the average distance between conjugate points provides a relevant correlation time scale. It can be 
proved that [58,64] if the sectional curvature K is bounded asO<L<K<H, then the distance d between two 
successive conjugate points is bounded by d > . The upper bound H of the curvature can then be approximated 
in our framework by 

~ fco + (Tfe , (69) 

so that we can define ri as (remember that dt = ds) 

n = d, = -j^=. (70) 

This time scale is expected to be the most relevant only as long as the curvature is positive and the fluctuations are 

small, compared to the average. 

Another time scale, referred to as T2, is related to the local curvature fluctuations. These will be felt on a length 
scale of the order of, at least, I = l/^/ok (the average fluctuation of curvature radius). The scale I is expected to 
be the relevant one when the fluctuations are of the same order of magnitude as the average curvature. Locally, the 
metric of a manifold can be approximated by [58] 

9ik - 6ik - ^Rikjiu''u'' , (71) 

where the are the components of the displacements from the point around which we are approximating the metric. 
When the sectional curvature is positive (resp. negative), lengths and time intervals — on a scale / — are enlarged 
(resp. shortened) by a factor {PK/6), so that the period has a fluctuation amplitude d2 given by d2 = ^^:^f=) 
replacing K by the most probable value ko one gets 

_ _l^ko 27r ^kl^' 
T2=d2 = — = ~ . (72) 

D ^ko (Tk 

Finally r in Eq. (68) is obtained by combining ti with T2 as follows 

r-i=rr^+r2-V (73) 

The present definition of t is obviously by no means a direct consequence of any theoretical result, but only a 
reasonable estimate. Such an estimate might well be improved independently of the general geometric framework. 

Now that all the quantities entering Eq. (65) have been fully defined, we can turn to item (iv), i.e., to the solution 
of Eq. (65). Whenever fc(f) has a non- vanishing stochastic component, any solution -(/"(i) has an exponentially growing 
envelope [65] whose growth-rate provides a measure of the degree of instability. How can one relate such a growth-rate 
with the Lyapunov exponent of the physical system? Let us recall that, for a standard Hamiltonian system of the 
form (1), the Lyapunov exponent can be computed as the following limit (see Appendix C, Eq. C12): 

A = lim 1 log erii) + -+e.ii)^m^-^^%ii) 

where the ^'s are the components of the tangent vector, i.e., of the pertiubation of a reference trajectory, which obey 
the tangent dynamics equation (25). In the case of Eisenhart metric, each component of the Jacobi vector field J can 
be identified with the corresponding component of the tangent vector ^; moreover, in Eq. (65) stands for any of the 
components of J, which obey the same effective equation. Thus, Eq. (74) becomes 

A=limllog^!l^l±im, (75) 
t^oo 2t ^ V^(0) + V2(0) 

where V'(i) is solution of Eq. (65). Equation (75) is our estimate for the (largest) Lyapunov exponent. 
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As a stochastic difFcrcntial equation, the solutions of Eq. (65) are properly defined after an averaging over the 
realizations of the stochastic process: referring to such an averaging as (•), we rewrite Eq. (75) as 

A^liml log + , (76) 

The evolution of (V'^), w^) ^-^d (i/'V'): i-e-, of the vector of the second moments of -0, obeys the following equation, 
which can be derived by means of a technique, developed by Van Kampen and sketched in Appendix D: 

d ( ^^^^ \ / 2 
It 




(V-^) = U^r -2fco 1 I (V-^) 1 (77) 



where fco and crfc are the mean and the variance of k{f), defined in Eqs. (66a) and (66b), respectively. Equation (77) 
can be solved by diagonalizing the matrix on the r.h.s. of (77). The result for the evolution of ('^^) + ('0^) is 

{i,\t)) + {i,\t)) = ((V'^(O)) + (V''(0))) exp(at) , (78) 

where a is the only real eigenvalue of the matrix. According to Eq. (76), the Lyapunov exponent is given by A = a/2, 
so that, by computing explicitly a, one then finds the final expression 

A(fco,<7fc,r) = i('A-^] , (79a) 




(79b) 



All the quantities fco, and T{ko,ak) can be computed as static averages, as functions of the energy per degree of 
freedom, e (see Eqs. (66a) and (66b)). Therefore — within the limits of validity of the assumptions made above — 
Eqs. (79) provide an approximate analytic formula to compute the largest Lyapunov exponent independently of the 
numerical integration of the dynamics and of the tangent dynamics. 

Let us remark that expanding Eqs. (79) in the limit (Tfe <^ ko one finds that 

Xoc al (80) 

which shows how close the relation is between curvature fluctuations and dynamical instability. 



B. Some applications 

Let us now discuss briefly the results of the application of the geometric techniques described up to this point to 
some Hamiltonian models. In particular, we shall consider two cases: a chain of coupled nonlinear oscillators (the 
so-called FPU (3 model, first introduced by Fermi, Pasta and Ulam in Ref. [66]) and a chain of coupled rotators (the 
1-d XY model). The reason of the choice of these two particular models is that they allow fully analytic calculations 
and are well-suited to show advantages and limitations of the theory. The geometric theory developed above has 
already been applied to many other cases, some of which will be addressed in Sec. V. For other applications we 
refer to the literature: in particular, a model of a homopolymer chain has been studied in Ref. [32], a model of a 
three-dimensional Lennard- Jones crystal has been studied in Ref. [34], and a classical lattice gauge theory has been 
considered in Ref. [39]. Without entering into the details, we would like to single out a result which is shared by all 
the models considered until now. In all these models the functional dependence of the largest Lyapunov exponent on 
the energy per degree of freedom s, in the low-sr limit, is numerically found to be 

A(£) cx . (81) 

No explanation of this "universal" behavior is yet at hand, and for some cases doubts about the validity of such a 
scaling with energy have been raised, because the numerical determination of Lyapunov exponents at low e is difficult. 
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However, the application of the geometric theory has provided a theoretical confirmation of this behavior in all the 

cases considered. 

The systems we now consider are 1-d models with nearest-neighbor interactions whose Hamiltonians H have the 
standard form (1) with 



N 



V = ^v{qi-q,-,) . (82) 
i=l 

The interaction potentials are, respectively, 

v{x) = ^x^ + ^x'^ (FPU-/3 model) , (83a) 

or 

t;(a;) = -Jcosa; (1-d XF model) . (83b) 

In Eq. (83a) we used u instead of the customary /3 in order to avoid confusion with the inverse temperature 0. We 

assume u > 0. 

The geometric quantities — in the framework of the Eisenhart metric — which are relevant in the quasi-isotropy 
approximation to describe the "effective" structure of the mechanical manifold, and which enter the geometric formula 
for the Lyapunov exponent, are the average and the root mean square (r.m.s.) fluctuations of the Ricci curvature of 
the mechanical manifold. They arc defined as statistical averages computed in the microcanonical ensemble (see Eqs. 
(66)). First, we will show how these microcanonical quantities can be computed starting from the canonical partition 
function, which can be calculated exactly for an infinite chain, i.e., N ^ oo, for both models (83). Then, we will 
apply this procedure to each of the two models (83) . 

The average and fiuctuations, within the microcanonical ensemble, of any observable function f{q), can be computed 
as follows, in terms of the corresponding quantities in the canonical ensemble. The canonical configurational partition 
function Z{(3) is given by 

Z{P) = J dqe-^^^'i^ (84) 

where dq = OiLi ^Qi- The canonical average (/)can of the observable / can be computed as 

(/)can = [Zm-' J dq f{q) e-''^(«) . (85) 

From this average, we can obtain the microcanonical average of /, (/)^, in the following (implicit) parametric form 
[67] 



(/)m(/3) = (/)can(/3) 



(/)^(^) (86) 



Note that Eq. (86) is strictly valid only in the thermodynamic limit; at finite A'', (/)^(/3) = (/)can(/3) + 

Contrary to the computation of (/) , which is insensitive to the choice of the probability measure in the N ^ oo 

limit, computing the fluctuations of /, i.e., of (5^/) = jfHf — {f))^), by means of the canonical or microcanonical 
ensembles yields different results. The relationship between the canonical — i.e. computed with the Gibbsian weight 
e"''^ — and the microcanonical fluctuations, is given by the Lebowitz-Percus-Verlet formula [67] 



(<5V)^(£) = (<5V)can(/3)-f 



2 ro, .X ,n^-,2 



5(/)can(/3) 



df3 



(87) 



where 



^ _ 0' g(^)can .... 
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is the specific heat at constant volume and (3 = P{e) is given in impHcit form by the second equation in (86). 

The average ko and the fluctuations ak of the Ricci curvature per degree of freedom are then obtained by replacing 
/ with the explicit expression for Ricci curvature, which, according to the definition given in Eq. (64), is 



N 



& 



(89) 



in Eqs. (86) and (87), respectively. 

We now turn to the two applications mentioned above. 



1. FPU 13 model 



For the FPU /3 model the dynamical observable which corresponds to the Ricci curvature reads, according to Eq. 
(89), 



N 



KR = 2N + 6u^{qi+i-qi) 



(90) 



Note that Kn is always positive and that this is also true for the sectional curvature along a physical geodesic. 
Computing the micro canonical average of Kr according to Eq. (86) we find that in the thermodynamic limit ko{s) is 
implicitly given by (the details are reported in Ref. [31]) 



g I 3 D-s/2ie) 
D_,/2{0) 



e{0) = 



8(7 



3 1 D-3/2iO) 



ko{e) , 



(91) 



where the Di, are parabolic cylinder functions [61] and ^ is a parameter proportional to /3, so that 6 G [0, +oo]. 
Let us now compute the fluctuations 



N 



(92) 



According to Eq. (87), first the canonical fluctuation, ((5^fcK)can(/3) = jf{{KR~ (-f^_R))^)can(/?), has to be computed 
and then a correction term must be added. For the canonical fluctuation we obtain [31] 



and the final result for the fluctuations of the Ricci curvature is 



D-3/2i9) 



(93) 



cv{0) 



dp 



2 \ 



6(6) 



1 



1 ^-3/2(0) 



(94) 



.02 D_y^{e)_ 

where (5^fcij)can(0) is given by (93), d{kR){0)/dP is given by 

djkRm ^ 3 0D'-s/2i^) + ^{0' - i)D-i/2{0)D-s/2{0) - '^OD^/M 

dp 8/103 ^-l/2(^) 

and the specific heat per particle cy is found to be [68,31] 



(95) 
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- 0275_3/2(0)[20D_i/2(0)+D_3/2W]} . 



(96) 



The microcanonical averages and fluctuations computed in Eqs. (91) and (94) are compared in Figs. 1 and 2 with 
their corresponding time averages computed along numerically simulated trajectories of the FPU /3-model with the 

potential (83a) for = 128 and = 512 with u = 0.1. Though the microcanonical averages have to be computed in 
the thermodynamic limit, the agreement between time and ensemble averages is excellent already at iV = 128. 




FIG. 1. Average Ricci curvature (Eisenhart metric) per degree of freedom, feo, vs. energy density e for the FPU-/? model. The 
continuous lino is the analytic computation according to Eq. (91); circles and squares are time averages obtained by numerical 
simulations with A'" = 128 and N = 512 respectively; u = 0.1. From Ref [24]. 




FIG. 2. Fluctuations of the Ricci curvature (Eisenhart metric), at vs. energy density e for the FPU-/3 model. Symbols and 
parameters as in Fig. 1; the continuous line now refers to Eq. (94). From Ref [24]. 

Before we comment on these results, we remark here that in many Hamiltonian dynamical systems different dynam- 
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ical regimes can be found as the energy per degree of freedom e is varied (sec, for a review, Ref. [69] and references 
quoted therein). In particular, in the FPU-/? model, a weakly chaotic regime is found for specific energies smaller than 
£c ~ 0.1/w, [70,71,24]. Although in the weakly chaotic regime the dynamics is chaotic (i.e., the Lyapunov exponent 
is positive, though small), mixing is very slow, as witnessed by the existence of a rather long memory of the initial 
conditions, i.e., of long relaxation times if the initial conditions are far from equilibrium. For e larger than Ec the 
dynamics is strongly chaotic and relaxations are fast. The precise origin of these phenomena is still to be understood. 
However, the geometric approach described here is able to provide a suggestive interpretation [24,33]. Let us consider 
Fig. 3, where the ratio of the fluctuations and the average curvature ak/ko is reported. As e ^ 0, crfc <C ko, so that the 
manifold looks essentially like a constant curvature manifold with only small curvature fluctuations. This situation 
corresponds to the weakly chaotic dynamical regime. On the contrary, as e is larger than Ec, ffe/Zeo tends to saturate 
towards a value of order unity, thus indicating that in the high-energy (strongly chaotic) regime the curvature fluctu- 
ations are of the same order of magnitude as the average curvature, so that the system no longer "feels" the isotropic 
(and integrable) limit. Hence the geometric approach can give a hint for understanding, at least qualitatively, the 
origin of weak and strong chaos in the Fermi-Pasta-Ulam model. 




10-2 10-1 101 10== 10= 



e 

FIG. 3. Fluctuations of the Ricci curvature (Eisenhart metric) divided by the average curvature, ak/ko, vs. energy density 
£ for the FPU-/9 model. 

The geometric theory also allows us to make a quantitative prediction for the Lyapunov exponent as a function 
of ko and ah via Eq. (79) , which turns out to be extremely accurate. The analytic result is shown in Fig. 4 and is 
compared with numerical simulations made for different values of A'', for the FPU-/? case in a wide range of energy 
densities — more than six orders of magnitude [25,31]. The agreement between theory and simulations is remarkably 
good. 
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FIG. 4. Lyapunov exponent A vs. energy density e for the FPU-/? model with u = 0.1. The continuous Une is the theoretical 
computation according to Eq. (79), while the circles and squares are the results of numerical simulations with N respectively 
equal to 256 and 2000. From Ref. [31]. 



2. 1-d XY model 



If the canonical coordinates qj and p,; arc given the meaning of angular coordinates and momenta, the 1-d XY 
model, whose potential energy is given in Eq. (83b), describes a linear chain of N rotators constrained to rotate on a 
plane and coupled by a nearest-neighbor interaction. This model can be formally obtained by restricting the classical 
Heisenberg model with 0(2) symmetry to one spatial dimension. The potential energy of the 0(2) Heisenberg model is 
V = —JJ2{i j) ' ^ii where the sum is extended only over nearest-neighbor pairs, J is the coupling constant and each 
s, has unit modulus and rotates in the plane. To each "spin" s, = {cos qt, sin qi), the velocity Sj = {— qi sin qi,qi cos qt) 
is associated, so that H = ^sf — JJ2{ij) ' ^j- 

This Hamiltonian system has two integrable limits. In the small energy limit it represents a chain of harmonic 
oscillators, as can be seen by expanding the potential energy in a power series, 



i=i 



(97) 



where Pi ^ q-i, whereas in the high-energy limit a system of freely rotating objects is found, because the kinetic energy 
becomes much larger than the bounded potential energy. 

The dynamics of this system has been extensively studied recently [68,72,28]. Numerical simulations and theoretical 
arguments independent of the geometric approach (see in particular Ref. [72]) have shown that also in this system 
there exist weakly and strongly chaotic dynamical regimes. It has been found that the dynamics is weakly chaotic in 
the low- and high-energy density regions, close to the two integrable limits. On the contrary, fully developed chaos is 
found in the intermediate-energy region. 

According to Eq. (89), the expression of the Ricci curvature Kr, computed with the Eisenhart metric, is 

N 

KR{q) = 2 J ^ cos(gi+i - qi). (98) 

We note that for this model a relation exists between the potential energy V and Ricci curvature Kr: 

V{q) = JN-^^^. (99) 
The average Ricci curvature can be again expressed by implicit formulas (see Ref. [31] for details) 
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ko{e) , 



(100) 



where the /,y's are modified Bessel functions of index [61]. The fluctuations are given by the impUcit equations 

4 J fiJI^ifSJ) - Io{f3J)h{(3J) - I3JIKI3J) 



liiPJ) l + 2{pjf -2pJh{pj)Io{pj)-2[l3Jh{(3J)f 



(101) 



In Figs. 5 and 6 a comparison between analytical and numerical results is provided for the average Ricci curvature 
and its fluctuations. The agreement between ensemble and time averages is again very good. 




FIG. 5. Average Ricci curvature (Eisenhart metric) per degree of freedom feo vs. specific energy e for the coupled rotators 
model. The continuous line is the result of an analytic computation according to Eq. (100); the full circles are time averages 
obtained by numerical simulations with N = 150; J = 1. Prom Ref. [31]. 
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FIG. 6. Fluctuation of the Ricci curvature (Eisenhart metric) vs. specific energy e for the coupled rotators model. Symbols 
and parameters as in Fig. 5; the continuous line now refers to Eq. (101). Prom Ref. [31]. 



Looking at Figs. 5 and 6 we realize that the low-energy weakly chaotic region has the same geometric properties as 
the corresponding region of the FPU model, as expected, since the two low-cncrgy intcgrablc limits arc the same. On 
the contrary, in the high-energy weakly chaotic region the fluctuations are far from being small with respect to the 
average curvature. The average curvature ko{e) vanishes as £ ^ oo. In this case the weakly chaotic dynamics seems 
related to the fact that the manifold (M x R^, g^;) looks almost flat along the physical geodesies. The boimds of the 
two weakly chaotic regions, as estimated in Ref. [72], coincide with the values of e where the asymptotic behaviours 
of k (low-energy region) and Uk (high-energy region) set in, respectively. Moreover, the case of the coupled rotators is 
very different from the FPU case, since the sectional curvature K[s) along a geodesic can take negative values. The 
probability P{e) that K{s) < can be analytically estimated in the following simple way. The explicit expression of 
the sectional curvature i^(7, ^), relative to the plane spanned by the velocity vector 7 = dq/dt and a generic vector 
^±7, is (see Appendix A, Eq. A43) 



OiO/c ■ 



d(f__^d(f_^ ^ d'^V 



so that, computing d'^V/dq^dq'' using the explicit form of V{q) given in Eq. (83b), we get 



J 



N 



^cos(g.+i-g.) [f+'-e]' 



(102) 



(103) 



for the 1-d XY model. We realize, by simple inspection of Eq. (103), that the probability of finding if < along 
a geodesic must be related to the probability of finding an angular difference larger than | between two nearest- 
neighboring rotators. Prom Eq. (103) we see that N orthogonal directions of the vector ^ exist such that the sectional 
curvatures relative to the TV planes spanned by these vectors together with 7 — are just cos((jj+i — qi). These 
directions are defined by the unit vectors of components (1,0,..., 0), (0, 1, 0, . . . , 0), . . . , (0, . . . , 0, 1). Hence the prob- 
ability P{e) of occurrence of a negative value of a cosine is used to estimate the probability of occurrence of negative 
sectional curvatures along the geodesies. This probability function, calculated using the Boltzmann weight, has the 
following simple expression [28,31] 



P{e) = 



/!^0(-cosa;)( 



,/3 J cos X 



dx // e^-^^°''^(ia; 



J—w 



277 Jo (/3 J) 



(104) 



where 0(.t) is the Heavyside unit step fimction and Iq the modified Bessel function of index 0. P{e) is plotted in Fig. 
7. We see that in the strongly chaotic region such a probability starts to increase rapidly from a very small value, 
while it approaches an asymptotic value P{e) ~ 1/2 when the system enters its high-energy weakly chaotic region. 
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FIG. 7. Estimate of the probability P(e) of occurrence of negative sectional curvatures in the 1-d XY model according to 
Eq. (104); J = 1. From Ref. [31]. 

When the sectional curvatures are positive'' chaos is produced by curvature fluctuations, hence we expect chaos to 

be weak as long as (7k /ka <C 1, and to become strong when (7^ w fco. On the contrary, when K{s) can assume both 
positive and negative values, the situation is much more complicated, for there are now two different and independent 
sources of chaos: negative curvature which directly induces a divergence of nearby geodesies, and the bumpiness of 
the ambient manifold which induces such a divergence via parametric instability. The results for the coupled rotators 
model suggest that as long as the negative curvatures are "few" they do not dramatically change the picture, and 
may strengthen the parametrically generated chaos, while when their occurrence is equally likely as the occurrence of 
positive curvatures, the two mechanisms of chaos seem to inhibit each other and chaos becomes weak^. 



'^The sectional curvature is always strictly positive in the FPU /3 model; in the 1-d XY model, in the low energy region, 
negative sectional curvatures can occur, but have a very small probability. 

®The fact that the two mechanisms, when comparable, can inhibit rather than strengthen each other can be considered a 
"proof" of the fact that their nature is intrinsically different. A similar situation is found also in some billiard systems, where 
there are two mechanisms which can produce chaos: (i) defocusing, due to positively curved boundaries, and {it) divergence of 
the trajectories due to scatterings with negatively curved boundaries [73]. 
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FIG. 8. Lyapunov exponent A vs. energy density e for the 1-d XY model with J = 1. The continuous Hue is the theoretical 

computation according to Eq. (79), while full circles, squares and triangles arc the results of numerical simulations with A'^, 
respectively, equal to 150, 1000, and 1500. The dotted line is the theoretical result where the value of k entering Eq. (79) has 
been corrected according to Eq. (106) with a = 150. Prom Ref. [31]. 



Such a qualitative picture is consistent with the result of the geometric computation of A for the coupled rotator 
model. The result of the application of Eq. (79) to this model is plotted in Fig. 8 (solid line). There is agreement 
between analytic and numeric values of the Lyapunov exponent only at low and high energy densities. Like in the 
FPU case, at low energy, in the quasi-harmonic limit, we find A(e) oc s^. At high energy A(e) cx e^^/^; here A(£) is a 
decreasing function of e because for e ^ oo the systems is intcgrable. 

However, in the intermediate energy range our theoretical prediction underestimates the actual degree of chaos of 
the dynamics. This energy range coincides with the region of fully developed (strong) chaos. According to the above 
discussion the origin of the underestimation can be found in the fact that the role of the negative curvatures, which 
appears to strengthen chaos in this energy range, is not correctly taken into account. The sectional curvature K{s), 
whose expression is given by Eq. (103), can take negative values with non- vanishing probability regardless of the value 
of e, whereas, as long as e < J. this possibility is lost in the K^placement of K by the Ricci curvature, due to the 
constraint (99), which implies that at each point of the manifold 

fcfl(e) > 2(J-£) . (105) 

Thus our approximation fails to account for the presence of negative sectional curvatures at values of s smaller than 
J. In Eq. (103) the cosines have different and variable weights, (^*^^ — ■f*)'^, which make it in principle possible to 
find somewhere along a geodesic a ii" < even with only one negative cosine. This is not the case for ka where all 
the cosines have the same weight. 

Let us now show how the theoretical results can be improved. Our strategy is to modify the model for K{s) in 
some ejfective way wliic;li takes into acx'ouiit the just mentioned difficulty of kii{s) to adequately model K{s). This 
will be achieved by suitably "renormalizing" ko or to obtain an "improved" gaussian process which can better 
model the behavior of the sectional curvature. Since our "bare" gaussian model underestimates negative sectional 
curvatures in the strongly chaotic region, the simplest way to renormalize the gaussian process is to shift the peak 
of the distribution V{Kr) toward the negative axis to make the average smaller. This can easily be done by the 
following rescaling of the average curvature ko: 

This correction has no influence either when P{£) ~ (below e ~ 0.2) or when P{e) ~ 1/2 (because in that case 

{kii{e)) 0). The simple correction (106) makes use of the information we have obtained analytically, i.e., of the P{e) 
given in Eq. (104), and is sufficient to obtain an excellent agreement of the theoretical prediction with the numerical 
data over the whole range of energies, as shown in Fig. 8. The parameter a in (106) is a free parameter, and its value 
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is determined so as to obtain the best agreement between numerical and theoretieal data. The resuh shown in Fig. 8 
(dotted Une) is obtained with a = 150, but also very different values of a, up to a ~ 1000, yield a good result, i.e., no 
particularly fine tuning of a is necessary to obtain a very good agreement between theory and numerical experiment. 

C. Some remarks 

Before moving to the second Part of the paper, let us now comment about some points of the material presented in 
the first Part of this Report. In particular, we would like to clarify the meaning of some of the approximations made 
and to draw the attention of the reader to some of the points which are still open. A better understanding of these 
points could lead, in our opinion, to a considerable improvement of the theory, which is still developing and can by 
no means be considered as a "closed" subject. 

What has been presented in this Section has a conceptual implication that goes far beyond the development of a 
method to analytically compute Lyapunov exponents. Rather, the strikingly good agreement between analytic and 
numeric Lyapunov exponents obtained at the price of a restriction of the domain of applicability of the analytic 
expression worked out for A - has three main implications: 

(i) the local geometry of mechanical manifolds contains all the relevant information about (in)stability of Hamilto- 
nian dynamics; 

(m) once a good model for the local source of instability of the dynamics is provided, then a statistical-mechanical- 
like treatment of the average degree of instability of the dynamics can be worked out, in the sense that we do 
not need a detailed knowledge of the dynamics but, by computing global geometric quantities, obtain a very 
good estimate of the average strength of chaos; 

{in) due to the variational formulation of newtonian dynamics, the Riemannian-geometric framework a-priori seems - 
and actually seemed in the past (as we have pointed out in Section II) - the natural framework to investigate the 
instability of Hamiltonian dynamics, however no evidence was available at all to confirm such an idea until the 
above mentioned development took place. It is now evident that the efforts to improve the theory by expanding 
its domain of applicability are worthwhile. 

We must warn the reader though against "blind" applications of formula (79), i.e. without any idea about the 
fulfilment, by the Hamiltonian model under investigation, of the conditions under which it has been derived. 

Prom a more technical point of view, one of the central results we have presented so far is the possibility of deriving, 
from the Jacobi equation, a scalar equation (Eq. (52)) describing the evolution of the Jacobi field J for a geodesic 
spread on a manifold. We would like to stress that such a result, though approximate, applies to a wide class of 
Hamiltonian systems. In fact, the only hypothesis needed to get such an equation is the quasi- isotropy hypothesis, 
i.e., the assumption that Riju ~ K,{s){gik9ji — dudjk)- Loosely speaking, such an assumption means that, locally, 
the manifold can be regarded as isotropic, i.e., there is a neighborhood of each point where the curvature can be 
considered constant. This does not imply at all that there are only small-scale fluctuations. There can be fluctuations 
of curvature on many scales, provided that they are finite and there is a cutoff at a certain point. The only case in 
which such an assumption will surely not hold is when there are curvature fluctuations over all scales. As will become 
clear in the following, this might happen when the manifold undergoes a topological change, and for "mechanical 
manifolds" this might happen at a phase transition. 

Other approximations come into play when one actually wants to model /C(s) along a geodesic with a stochastic 
process. It is true that replacing the sectional curvature by the Ricci curvature requires that the fluctuations are not 
only finite, but also small. Moreover, we use global averages to define the stochastic process, and here it is crucial 
that the fluctuations do not extend over too large scales. Thus Eq. (58) has a less general validity than Eq. (52). A 
way to improve the theory might be to try to replace the sectional curvature with some quantity related also to the 
gradient of the Ricci curvature, in order to make the replacement of sectional curvature less sensitive to the large 
scale variations of the Ricci curvature. 

To get an explicit solution of Eq. (52), an even less general situation must be considered, through the following 
steps: 

(i) using the Eisenhart metric; 

(m) considering standard systems where the kinetic energy does not depend on the (/'s; 
{Hi) estimating the characteristic correlation time r of the curvature fluctuations. 
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As to item (Hi), wc have already remarked that our estimate given in Eq. (73) is by no means a consequence of any 
theoretical result, but only a reasonable estimate which could surely be improved. 

As to item (m), the case of a more general kinetic energy matrix 7^ 5ij, though not conceptually different, is 
indeed different in practice and the final result is not expected to hold in the same form for that case. 

Finally, (item (z)) should not reduce significantly the generality of the result. In fact, considering the Eisenhart 
metric only makes the calculations feasible, and in principle nothing should change, if one were able to solve Eq. 
(58) in the case of the Jacobi metric (see the discussion in Ref. [35]). However, Eisenhart and Jacobi metrics are 
equivalent for what concerns the computation of the average instability of the dynamics [35], but they might not 
be equivalent for other developments of the theory. This in view of the fact that (ME,gj) is a manifold which has 
better mathematical properties than (M x R^,^^;): {ME,gj) is a proper Riemannian manifold, it is compact, all of 
its geodesies are in one-to-one correspondence with mechanical trajectories, its scalar curvature does not identically 
vanish as is the case of (M x R^, g^;), it can be naturally lifted to the tangent bundle where the associated geodesic 
flow on the submanifolds of constant energy coincides with the phase space trajectories. 

Let us finally add a comment on the application of the theory to the calculation of the Lyapunov exponent for 
the one-dimensional XY model. We have seen that although the predictions of the theory compare reasonably well 
with the numerical simulations, there is an intermediate energy range in which a correction must be added. As will 
become clear in the second Part, the very first assumption (quasi-isotropy) should not be satisfied for this model, due 
to the presence of topology changes in the mechanical manifolds, in fact the difficulties of the theory begin just at the 
energy density which corresponds to the appearance of a large number of critical points of the potential energy (see 
next Sections). 

V. GEOMETRY AND PHASE TRANSITIONS 

In the previous Sections we have shown how simple concepts belonging to classical differential geometry can be 

successfully used as tools to build a geometric theory of chaotic Hamiltonian dynamics. Such a theory is able to 
describe the instability of the dynamics in classical systems consisting of a large number N of mutually interacting 
particles, by relating these properties to the average and the fluctuations of the curvature of the configuration space. 
Such a relation is made quantitative through Eq. (79), which provides an approximate analytical estimate of the 
largest Lyapunov exponent in terms of the above-mentioned geometric quantities, and which compares very well with 
the outcome of numerical simulations in a number of cases, two of which have been discussed in detail at the end of 
Sec. IV. 

The macroscopic properties of large- A'' Hamiltonian systems can be understood by means of the traditional meth- 
ods of statistical mechanics. One of the most striking phenomena that may happen in such systems is that when 
the external parameters (e.g., either the temperature or the energy) are varied until some critical value is reached, 
the macroscopic thermodynamical quantities may suddenly and even discontinuously change, so that, though the 
microscopic interactions are the same above and below the critical value of the parameters, its macroscopic prop- 
erties may be completely different. Such phenomena are referred to as phase transitions. In statistical mechanics, 
phase transitions are explained as true mathematical singularities that occur in the thermodynamic functions in the 
limit N ^ 00, the so-called thermodynamic limit^ [75]. Such singularities come from the fact that the equilibrium 
probability distribution in configuration space, which in the canonical ensemble is the Boltzmann weight 

&an(gi, . . . , gjv) = exp [-(3V{qi, gjv)] , (107) 

where (3 = l/ksT, V is the potential energy, and Z = f dqe~^^^'^^ is the configurational partition function, can itself 
develop singularities in the thermodynamic limit. 

The statistical-mechanical theory of phase transitions is one of the most elaborate and successful physical theories 
now at hand, and at least as continuous phase transitions are concerned, also quantitative predictions can be made, 
with the aid of renormalization-group techniques, which are in very good agreement with laboratory experiments and 
numerical simulations. We are not going to discuss this here, referring the reader to the vast literature on the subject 
[76-82]. 



® According to G. E. Uhlonbock [74], the use of the thermodynamic limit as an explanation of the singularities of the partition 
function was suggested for the first time by Kramers in the 1938 Leiden conference on Statistical Mechanics. 
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However, the origin of the possibiUty of describing Hamihonian systems via equihbrium statistical mechanics are the 
chaotic properties underlying the dynamics. In fact, though it is not possible to prove that generic Hamiltonian systems 
of the form (1) are ergodic and mixing, the fact that the trajectories arc mostly chaotic (i.e., for the overwhelmingly 
majority of the trajectories positive Lyapunov exponents are found) means that such systems can be considered 
ergodic and mixing for all practical purposes^°. 

The observation that chaos is at the origin of the statistical behaviour of Hamiltonian systems and that chaotic 
dynamics can be described by means of the geometric methods described above, naturally leads to the foUwing two 
questions: 

1. Is there any specific behaviour of the Lyapunov exponent^^ when the system undergoes a phase transition? 

2. What are the geometric properties of the configuration space manifold in the presence of a phase transition? 

The aim of the present Section is to discuss these two questions. We shall first give a phenomenological description 
which follows from numerical experiments, and then we shall concentrate on the particular case of the mean-field XY 
model where the geometrical quantities can be analytically calculated. From the discussion of these questions and 
from the (at least partial) answers that we find, we are lead to put forward a Topological Hypothesis about phase 
transitions, which will be discussed in Sec. VI. 



In order to look for an answer to question 1 above, we now review the numerical results that have been obtained 
until now for various Hamiltonian dynamical systems which show a phase transition when considered as statistical- 
mechanical models for macroscopic systems in thermal equilibrium. 

The first attempt to look for a chaotic-dynamic counterpart of an equilibrium phase transition is in the work by 
Butcra and Caravati (BC) [83]. BC considered a two-dimensional XY model, i.e., a Hamiltonian dynamical system 
of the form (1) with the potential energy 



where the i^'s are angles, i and j label the sites of a square lattice and the sum runs over all the nearest-neighbor sites. 
This model is the two-dimensional version of the one studied in Sec. IV B 2. As the temperature is decreased, such a 
system undergoes a peculiar phase transition (referred to as the Bcrczinskij-Kosterlitz-Thouless, or BKT, transition) 
from a disordered phase to a quasi-ordered phase where, though no true long-range order is present, correlation 
functions decay as power laws, as occurs at a critical point [82]. Since there are no singularities in the finite-order 
derivatives of the free energy, the BKT transition is sometimes classified as an "infinite-order" phase transition. BC 
computed the Lyapunov exponent A as a function of the temperature, and found that A(T) followed a rather smooth 
pattern; however, in a region around the transition, the dependence of A on T changed from a steeply increasing 
function to a much less steep one. 

BC's pioneering paper has been the only example of a study of this kind for a long period. However, very recently 
there has been a renewed interest in the study of the behaviour of Lyapunov exponents in systems undergoing phase 
transitions, and a number of papers appeared [36-38,84-94]. The two-dimensional XY model has been reconsidered, 
together with the three-dimensional case, in Ref. [36]. We remark that in three spatial dimensions the XY model 
undergoes a standard continuous (second-order) phase transition accompanied by the breaking of the 0{2) symmetry 
of the potential energy (108). The behaviour of the Lyapunov exponent A as a function of the temperature T is shown 
in Figs. 9 and 10. 



^°This has boon rocontly cxtondod to non-equilibrium statistical mechanics via the "Chaotic Hypothesis" [62]. 
'^'^As in the first part of the paper, we consider only the largest Lyapunov exponent, which is referred to as just the Lyapunov 
exponent. 



A. Chaotic dynamics and phcise transitions 
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FIG. 9. Lyapunov exponent A vs. the temperature T for the two-dimensional XY model: the circles refer to a 10 x 10, the 
squares to a 40 x 40, the triangles to a 50 x 50, and the stars to a 100 x 100 lattice, respectively. The critical temperature of 
the BKT transition is Tc — 0.95 and is marked by a dotted vertical line. From Ref. [36]. 




FIG. 10. Lyapunov exponent A vs. the temperature T for the three-dimensional XY model, numerically computed on a 
A'' = 10 X 10 X 10 lattice (solid circles) and on a A'^ = 15 x 15 x 15 lattice (solid squares), respectively. The critical temperature 
of the phase transition is Tc ^ 2.15 and is marked by a dotted vertical line. From Ref. [36]. 

The behavior found for the two-dimensional case confirms the BC results. The three-dimensional case shows a 
similar behavior, but the change of the shape of the A(T) function near Tc is somehow sharper than in the previous 
case. 

Dellago and Posch [88] considered an extended XY model, whose potential energy is 



y = 2 - 2 ^ cos 

{iJ) 



ipi - 



(109) 



which includes the standard XY model (108) for = 2. On a two-dimensional lattice the transition, which is a 
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continuous BKT transition for = 2, becomes a discontinuous transition as ~ 100. The results for the Lyapunov 
exponent A show that for any considered vahie of there is a change in the shape of A(T) close to the critical 
temperature. 

One of the systems which have received considerable attention in this framework is the so-called lattice if'^ model, 
i.e., a system with a Hamiltonian of the form (1) and a potential energy given by 



V 



J 



ihj) 



(110) 



where the (p's are scalar variables, (fi G [— oo,+oo], defined on the sites of a d-dimensional lattice, and and u are 
positive parameters. The lattice </j^ model has a phase transition at a finite temperature provided that d > 1. The 
existence of such a transition, which belongs to the universality class of the d-dimensional Ising model, can be proved 
by renormalization-group arguments (see e.g. [77,95]). The cases d = 2 and d = 3 have been considered in Refs. [38] 
and [37], respectively. Moreover, in Ref. [37] also some vector versions of this model have been considered, namely, 
systems with potential energy given by 



4! 



H 2 



E(<)^ 



(111) 



where a runs from 1 to n, labelling the components of the vectors ipi = [ipj, . . . , (/j"). The potential energy (111) is 
0(n)-invariant; in the case n = 1 we recover the scalar model (110). Figures 11 and 12 show the behaviour of A in 
the (^^ model, in two and three dimensions, respectively. 
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FIG. 11. Lyapunov exponent A vs. the energy per particle e, numerically computed for the two-dimension al 0(1) model, 
with N = 100 (solid circles), A^" = 400 (open circles). A'' = 900 (solid triangles), and A'' = 2500 (open triangles), respectively. 
The critical energy is marked by a vertical dotted line, and the dashed line is the power law e^. Prom Ref. [38]. 
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FIG. 12. Lyapunov exponent A vs. the temperature T for the three-dimensional ip^ model. Pull circles correspond to the 
0(1) (scalar) case, open circles to the 0(2) case, and open triangles to the 0(4) case. Prom Ref. [37]. 



Again we see that the Lyapunov exponent is sensitive to the presence of the transition, and that the shape of A(T) 
close to the transition is highly model-dcpcndcnt. Moreover, such a shape can be significantly different within the 
same model as its parameters are varied. For instance, in the tp'^ model, A can be either a monotonously increasing 
function of T or can display a maximum close to Tc, depending on the values of and u [38]. 

The Lyapunov exponents of systems undergoing phase transitions of the solid-liquid type have been recently deter- 
mined numerically: Dellago and Posch (DP) considered, in two dimensions, a system of hard disks [86], a Lorcntz- 
gas-like model and a Lennard- Jones fluid [87], and, in three dimensions, a system of hard spheres [89]. DP found that 
in all these systems the Lyapunov exponent is sensitive to the phase transition, and again the shape of A is different 
for different models, the common feature being that A attains a maximum close, if not at, the transition. Similar 
results have been obtained by Mehra and Ramaswamy [90] . Bonasera et al. [84] considered a classical model of an 
atomic cluster, whose particles interact via phenomenological pair potentials of the form 

t;(r) =ae-(^) -c(^)'^ , (112) 

and of a nuclear cluster, with nucleons interacting via Yukawa pair potentials. Such systems undergo a so-called 
"multifragmentation" transition at a critical (model-dependent) temperature Tc- Bonasera et al. computed the 
Lyapunov exponents of these systems by means of numerical simulations at different temperatures. The resulting 
A(T) of both systems develops a sharp maximum close to T^. 

The numerical evidence that we have reviewed above clearly shows that the Lyapunov exponent of a Hamiltonian 
dynamical system is sensitive to the presence of a phase transition. However, the interpretation of the observed 
behavior as it now stands is very difficult, because each model behaves differently and the behavior of A close to 
the transition docs not apparently exhibit any universal feature: on the contrary, the shape of A(r) can depend on 
the values of the parameters of the model. Moreover, the qualitative behavior of A(r) appears to be only weakly 
dependent on whether the transition is accompanied by a symmetry breaking or not, as in the case of the XY model: 
the shape of A(T) in two dimensions, where there is not any breaking of the 0(2) symmetry of the potential energy 
below the BKT transition temperature, is not dramatically different from that of the three-dimensional case where 
the phase transition is accompanied by a symmetry breaking. In the latter case the "knee" of the A(r) curve is 
sharper, but it would be difficult to discriminate between the two cases only by looking at the A(T) curve. Therefore, 
though clearly sensitive to the presence of a phase transition, the Lyapunov exponent does not seem a "good" probing 
observable for the occurrence of a symmetry-breaking phase transition. 
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B. Curvature and phase transitions 



In Section IV we have seen that the origin of chaos in Hamiltonian mechanics can be understood from a geometrical 
point of view, and that the Lyapunov exponents arc closely related to a geometric quantity, i.e., to the fluctuations of 
the Ricci curvature of the configuration space. Thus, it is natural to investigate whether such a geometric observable 
also has some peculiar behaviour close to the phase transition. As we shall see in the following, the fluctuations of the 
curvature do indeed have such a peculiar behaviour which, in turn, suggests a topological intepretation of the phase 
transition itself. 

The Ricci curvature along a geodesic of the enlarged configuration spacetime equipped with the Eisenhart metric, 
which we denoted by Kr in the previous Sections, is given by the Laplacian of the potential energy — see Eq. 64. In 
the case of the XY model we obtain, as already shown in §IVB 2, 

Kr = 2N -2V = 2^cos{^pi- ^pj) . (113) 

The root mean square fluctuation of Kr divided by the number of degrees of freedom N ^ i.e., 

<yk=(^^{Kl)-{Knf^"' , (114) 
is plotted in Figs. 13 and 14 for the 2-d and 3-d cases, respectively. 
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FIG. 13. Fluctuations of the Ricci curvature (Eisenhart metric), ak{T), vs. the temperature T for the two-dimensional XY 
model. The solid circles axe numerical values obtained for a 40 x 40 lattice; the dashed line is only a guide to the eye. The 
critical temperature of the BKT transition is Tc ^ 0.95 and is marked by a dotted line. From Ref. [36]. 
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FIG. 14. As Fig. 13, for the three-dimensional XY model. Here AT = 10 x 10 x 10, and the critical temperature of the phase 
transition is Tc ~ 2.15. From Ref. [36]. 

In the case of the (p^ model with 0{n) symmetry, the Ricci curvature Kji is given by [37,38] 



n N 



n N 



Nn{2Jd - r') + A(n + 2) ^ Y.{^?f • 



(115) 



a — 1 i—1 



The r.m.s. fluctuation of Kr, at, is plotted against the energy per degree of freedom, e, in the case of the two- 
dimensional 0(1) ip^ model in Fig. 15, and against the temperature T in the case of the two-dimensional 0{2) ip"^ 
model in Fig. 16, and for the three-dimensional 0{n) ip^ models in Fig. 17. 
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FIG. 16. Curvature fluctuations ak vs. the temperature T for the two-dimensional 0(2) (p'^ model, numerically computed 
on a square lattice of 30 x 30 sites. The critical temperature Tc of the BKT transition is located at Tc ~ 1.5. Prom Ref. [37]. 
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FIG. 17. Curvature fluctuations ak vs. the temperature T for the three-dimensional ip"^ model. Full circles correspond to the 
0(1) (scalar) case, open circles to the 0(2) case, and open triangles to the 0(4) case. From numerical simulations performed 
on an 8 X 8 X 8 cubic lattice, reported in Ref. [37]. 

Looking at Figs. 13-17, one can clearly see that when a symmetry-breaking phase transition occurs, a cusp-like 
("singular") behavior of the curvature fluctuations is found at the phase transition point (Figs. 14, 15 and 17), while, 
when only a BKT transition is present, no cusp- like pattern is observed^^ (Figs. 13 and 16). We can summarize 
these results by saying that, in general, curvature fluctuations always show a cusp-like behavior when a continu- 



Although the cusp-like behaviour is lost, indeed some change of behavior is still visible in Figs. 13 and 16 close to the critical 
temperature, so that a BKT transition appears as "intermediate" between the absence of a phase transition and the presence 
of a symmetry-breaking phase transition. 
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ous symmetry-breaking phase transition is present, and, within numerical accuracy, the cusp occurs at the critical 
temperature. No counter-examples have yet been found to this general rule. 

The fact that the Lyapunov exponent is sensitive to the phase transition can now be understood, in the light of the 
fact that, as shown in the first part of the present Report, chaos can be described geometrically and, under rather 
general assumptions, the Lyapunov exponent is closely related to the fluctuations of the Ricci curvature (see Eq. 80). 
However, contrary to the Lyapunov exponent, which, although sensitive to the phase transition, is not a good probing 
observable of the presence of a symmetry-breaking phase transition, the fluctuation of the Ricci curvature, (ik, is 
a good probing observable, for, as plotted as a function of the temperature, it exhibits a clearly peculiar ("cuspy") 
pattern when a symmetry-breaking phase transition is present, and a rather smooth pattern otherwise. The difference 
among A(T) and <Jk{T) as probing observables of the phase transition can be appreciated by comparing Figs. 12 and 
17, where A(T) and (ti~{T) are reported, respectively, for the 0{n) models. In Fig. 12, the A(T) curves for different 
values of n are qualitatively different, while in Fig. 17 the crfe(T) curves look strikingly similar, while being clearly 
different from the curve for the 2-d 0{2) case (Fig. 16), where only a BKT transition is present. The same can be 
said in the case of the 2-d and the 3-d 



At this point, it is worthwile to point out that we can apply the geometric formula (79) for the Lyapunov exponent to 
estimate A for all these models, since both fco and <Tk have been numerically computed. As shown in Refs. [36 38], one 
finds that in general, although the qualitative behavior of the Lyapunov exponent is well reproduced, the quantitative 
agreement between the values of A extracted from the numerical simulations and those obtained applying Eq. (79) is 
not good, in a neighborhood of the phase transition. 

However, this is to be expected, because among the assumptions under which the formula (79) was derived there 
was the hypothesis that the fluctuations of the curvature should be not too large, and this is clearly not true close to 
a phase transition, as we have just shown^^. 



The mean-fleld XY model [96] describes a system of N equally coupled planar classical rotators. It is defined by a 
Hamiltonian of the class (1) where the potential energy is 



Here (pi € [0, 2tt] is the rotation angle of the i-th rotator and h is an external field. Defining at each site i a classical 
spin vector Si = {cos ipi, sin ipi) the model describes a planar (XY) Heisenberg system with interactions of equal 
strength among all the spins. We consider only the ferromagnetic case J > 0; for the sake of simplicity, we set J = 1. 
The equilibrium statistical mechanics of this system is exactly described, in the thermodynamic limit, by mean-field 
theory [96]. In the limit ft. ^ 0, the system has a continuous phase transition, with classical critical exponents, at 
Tc = 1/2, or Ec = 3/4, where e = E/N is the energy per particle. 

The Lyapunov exponent A of this system is extremely sensitive to the phase transition. In fact, according to 
numerical simulations reported in Refs. [91,93,97,98], A(£) is positive for < e < Sc, shows a sharp maximum 
immediately below the critical energy, and drops to zero at £c in the thermodynamic limit, where it remains zero in 
the whole region e > £c, which corresponds to the thermodynamic disordered phase. In fact in this phase the system 
is integrable, reducing to an assembly of uncoupled rotators. These results are valid in the thermodynamic limit 
AT —> 00 in the sense that they have been obtained by estimating the infinite N limit of finite N numerical simulations 
[91,93]: in the whole region e > £c the Lyapunov exponent, numerically computed for systems with different numbers 
of particles A^, behaves as A oc A^~^/^, so that it extrapolates to zero at A^ ^ oo. 

These results have received a theoretical confirmation in recent work by M.-C. Firpo [92] based on the application 
of the geometric techniques described in the first part of the present paper. Firpo has computed analytically {ka) 



^^The results of the formula (79) can be improved using procedures which are specific of the model under consideration and 
which we are not going to describe here (see Ref. [36] for the XY case and Ref. [37] for the tp'^ case, respectively). 



1. Geometric estimate of the Lyapunov exponent 



C. The mean-field XY model 




(116) 
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and in the thermodynamic limit for the mean- field XY model, showing that such quantities indeed have a 

singular behavior at Sc (sec Fig. 18). Using these quantities and Eq. (79), Firpo has obtained the analytical estimate 
for A(£) reported in Fig. 19; it is remarkable that also the behavior A oc iV^^/"^ at e > £c has been extracted from 
this theoretical calculation. This result gives a theoretical confirmation to the qualitative behavior of the Lyapunov 
exponent extrapolated from the numarical simulations. Moreover, Firpo's analytical results are in good quantitative 
agreement with numerical results reported in Refs. [93,97], also close to the phase transition, at variance with the 
cases of the nearest-neighbor XY and ip^ models considered earlier. A tentative explanation of that the application 
of the geometric formula (79) gives such good quantitative results in the present case can be that the mean-field 
character of the model prevents the curvature fluctuations from being too wild. 
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FIG. 18. Mean field XY model: analytic expression for the microcanonical averages of the Ricci curvature (solid curve) and 
of its fluctuations (dot-dashed curve) . Prom Ref . [92] . 
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FIG. 19. Moan field XY model: analytic expression for the Lyapunov exponent (solid curve). The curves above the transition 
are finite- results for Af = 80 and = 200: here A oc N''^^^ . From Ref. [92]. 
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VI. PHASE TRANSITIONS AND TOPOLOGY 



In the previous Section we have reported results of numerical simulations for the fluctuations of observables of a 
geometric nature (e.g., configuration-space curvature fluctuations) related to the Ricmannian gcomctrization of the 
dynamics in configuration space^^. These quantities have been computed, using time averages, for many different 
models undergoing continuous phase transitions, namely (^^ lattice models with discrete and continuous symmetries 
and XY models. In particylar, when plotted as a function of either the temperature or the energy, the fluctuations of 
the curvature have an apparently singular behavior at the transition point. Moreover, we have seen that the presence 
of a singularity in the statistical-mechanical fluctuations of the curvature at the transition point has been proved 
analytically for the mean-field XY model. 

The aim of the present Section is to try to understand on a deeper level the origin of this peculiar behaviour. In 
§ VI A, wc will show, rising abstract geometric models, that a singular behaviour in the fluctuations of the curvature 
of a Riemannian manifold can be associated with a change in the topology of the manifold itself. By "topology 
change" we mean the following. Let us consider a surface Sg which depends on a parameter e in such a way that, 
upon varying the parameter, the surface is continuously deformed: as long as the different deformed surfaces can be 
mapped smoothly one onto another^^, the topology does not change; on the contrary, the topology changes if there 
is a critical value of the parameter, say Sc, such that the surface »Se>£^ can not any more be mapped smoothly onto 

The observation that a singularity in the curvature fluctuations of a Riemannian manifold, of the same type of those 
observed numerically at phase transitions, can be associated with a change in the topology of the manifold, leads us 
to conjecture that it is just this mechanism that could be at the basis of thermodynamic phase transitions. Such a 
conjecture was originally put forward in Ref. [36] as follows: a thermodynamic transition might be related to a change 
in the topology of the configuration space, and the observed singularities in the statistical-mechanical equilibrium 
measure and in the thermodynamic observables at the phase transition might be interpreted as a "shadow" of this 
major topological change that happens at a more basic level. We will refer to this conjecture as the Topological 
Hypothesis (TH). 

The remaining part of the present Section is devoted to a discussion of the TH and of its validity. In § VI B 1 we will 
report on a purely geometric, and thus still indirect, further indication that the topology of the configuration space 
might change at the phase transition, which has been obtained from numerical calculations for the model on a 
two-dimensional lattice [99,100]. Then, in §VIB 3, we will describe a direct confirmation of the TH, i.e., we will show 
that a topological change in configuration space, which can be related with a phase transition, indeed occurs in the 
particular case of the mean- field XY model [101]. Finally, in §VIC we will reformulate the TH in a more precise way, 
taking advantage of the previously discussed examples, and §VID will be devoted to a general discussion of the many 
points that are still open and on the future perspectives of the geometrical and topological approach to statistical 
mechanics. 



A. From geometry to topology: abstract geometric models 

Let us now describe how a singular behavior of the curvature fluctuations of a manifold can be put in correspondence 
with a change in the topology of the manifold itself. For the sake of clarity, we shall first discuss a simple example 
concerning two-dimensional surfaces [36,37], and then we will generalize it to the case of A''-dimensional hypersurfaces 
[99,100]. 

The simple geometric model wc are going to describe concerns surfaces of revolution. A surface of revolution S € 
is obtained by revolving the graph of a function / around one of the axes of a Cartesian plane, and can be defined, 
in parametric form, as follows [102] : 

S{u,v) = {x{u,v),y{u,v), z{u,v)) = {a{u) cosv, a{u) sin v,b{u)) , (H^) 

where either a{u) = f{u) and b{u) = u, if the graph of / is revolved around the vertical axis, or a{u) = u and 
b{u) = f{u), if the graph is revolved around the horizontal axis; in both cases, u and v are local coordinates on the 
surface S: v £ [0, 27r] and u belongs to the domain of definition of the function /. 



More precisely, we considered the enlarged configuration space-time, endowed with tlie Eisenhart metric. 
The different surfaces are then said to be diffeomorphic to each other (see Appendix B). 
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Let us consider now in particular the two families of surfaces of revolution defined as: 

J^s = if e{u) cosv, fs{u) smv,u) , 

and 

Qs = {u cosv,usmv, fs{u)) , 

where 



and £min = — J- Some cases are shown in Fig. 20. 



(118a) 
(118b) 
(119) 










FIG. 20. Some representatives of the two families of surfaces and Ge defined in Eqs. (118a) and (118b) respectively. Each 

family is divided into two subfamilies by the critical surface corresponding to Ec = (middle members in the picture) . Members 
of the same subfamily are diffeomorphic, whereas the two subfamilies are not diffeomorphic to each other. From Ref. [37]. 

There exists for both families of surfaces a critical value of £, Cc = 0, corresponding to a change in the topology of 
the surfaces: the manifolds are diffeomorphic to a torus for £ < and to a sphere for £ > 0; the manifolds 
Ge are diffeomorphic to t,wo spheres for £ < and to one sphere for £ > 0. Tlic Eulcr-Poincarc characteristic (see 
Appendix B, Eq. BIO) is xi^e) = if £ < 0, and xi^s) = 2 otherwise, while xiGe) is 4 or 2 for e negative or positive, 
respectively. 

Wc now turn to the definition and the calculation of the curvature fluctuations on these surfaces. Let M belong to 
one of the two families; its gaussian curvature K is [102] 



K 



a'{a"b' - b'a") 



l(&'2 



-J2\2 



(120) 



where a{u) and b{u) are the coefficients of Eq. (117), and primes denotes differentiation, with respect to u. The 
fluctuations of K can be then deflned as 



'K 



Jm 



/ KdS 

M 



(121) 



where A is the area of M and dS is the invariant surface element. Both families of surfaces exhibit a singular behavior 
in aK as £ ^ £c, as shown in Fig. 21, in spite of their different curvature properties on the average^^. 



^For instance, {K){£) = for J^g as e < 0, while for Gs the same average curvature is positive and diverges as e —^ 0. 
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FIG. 21. The fluctuation ax of the gaussian curvature of the surfaces J^e and Ge is plotted vs e. a is defined in Eq. (121), e 

is shifted by Emin = 0.25 for reasons of clarity of presentation, (a) refers to and (b) refers to J-g- The cusps appear at £ = 
where the topological transition takes place for both JT^ and Qe- From Ref. [37]. 

We are now going to show that the result we have just obtained for two-dimensional surfaces has a much more 
general validity: a generic topology change in an n-dimensional manifold is accompanied by a singularity in its 
curvature fluc;tuations [99]. In order to do that, we have to make use of some concepts belonging to Morse theory, 
which will also be used in §V1B3 below; the main concepts of Morse theory are sketched in Appendix B, where also 
references to the literature are given. 

We consider then a hypersurfacc of which is the m- level set of a function / defined in R-'^, i.e., a submanifold 
of R''^ of dimension n = N — 1 defined by the equation 



fixi, 



,Xn) = u 



(122) 



such a hypersurfacc can then be referred to as f^^{u). Let us now assume^''' that / is a Morse functAon, i.e., such that 
its critical points (i.e., the points of R^ where the differential df vanishes) are isolated. One of the most important 
results of Morse theory is that the topology of the hypersurfaces f~^{u) can change only crossing a critical level 
f~^{uc), i.e., a level set containing at least one critical point of /. This means that a generic change in the topology 
of the hypersurfaces can be associated with critical points of /. Now, the hypersurfaces f^^{u) can be given a 
Riemannian metric in a standard way [103] , and it is possible to analyze the behavior of the curvature fluctuations in 
a neighborhood of a critical point. Let us assume, for the sake of simplicity, that this critical point is located at = 
and belongs to the level Uc = 0. Any Morse function can be parametrized, in the neighborhood of a a;o, by means 



'^This is not a strong assumption: in fact, it can be shown that Morse functions are generic (see Appendix B). 
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of the so-called Morse chart, i.e., a system of local coordinates {yi} such that f{y) = /(.To) — X^iLi Vi + IZiLk+i vl 
{k is the Morse index of the critical point). Then standard formulae for the Gauss curvature K of hypersurfaces of 
R''^ [103] can be used to compute explicitly the fluctuations of the curvature, ck, of the level set f~^{u). Numerical 
results for the curvature fluctuations are reported in Fig. 22 and show that also at high dimension cr|. develops a 
sharp, singular peak as the critical surface is approached (for computational details, see Ref. [100]). 




FIG. 22. Fluctuations of the Gauss curvature of a hypersurface f^^iu) of R*'^ vs. u close to a critical point, o"^^ is 
reported because it has the same dimensions of the scalar curvature. Here dim(/~'^(w)) = 100, and the Morse indexes are: 
fe = 1, 15, 33, 48, represented by solid, dotted, dashed, long-dashed lines respectively. 



B. Topology changes in configuration space and phase transitions 

As we have discussed in Section V, the curvature fluctuations of the configuration space exhibit cusp-like patterns 
in presence of a second order phase transition. A truly cuspy pattern, i.e., an analytic discontinuity, is mathematically 
proven in the case of mean- field XY model. In §VIA, we have shown that singular patterns in the fluctuations of 
the curvature of a Riemannian manifold can be seen as consequences of the presence of a topology change. Hence, 
we are led to the Topological Hypothesis (TH), i.e., to conjecture that at least continuous, symmetry-breaking phase 
transitions are associated with topology changes in the configuration space of the system. 

However, an important question arises, in that the fluctuations of the curvature considered in Sec. V havc^ b(^en 
obtained as time averages, computed along the dynamical trajectories of the Hamiltonian systems under investigation 
(or as statistical averages computed analytically, as in the case of the mean-field XY model). Now, time averages of 
geometric observables are usually found to be in excellent agreement with ensemble averages [24,31,36-38] so that one 
could argue that the above mentioned singular-like patterns of the fluctuations of geometric observables are simply the 
precursors of truly singular patterns due to the fact that the measures of all the statistical ensembles tend to become 
singular in the limit N ^ oo when a phase transition is present. In other words, geometric observables, like any other 
"honest" observable, already at flnite A'' would feel the eventually singular character of the statistical measures, i.e., 
of the probability distribution functions of the statistical-mechanical ensembles. If this were the correct explanation, 
we could not attribute the cusp-like patterns of the curvature fluctuations to any special geometric features of conflg- 
uration space, and the cusp-like patterns observed in the numerical simulations could not be considered as (indirect) 
conflrmations of the TH. 

In order to elucidate this important point, three different paths have been followed: (i) purely geometric information 
about certain submanifolds of configuration space has been worked out independently of the statistical measiircs in 
the case of the two-dimensional </j^ model, and the results lend indirect support to the TH [99]; (ii) a direct numerical 
confirmation of the TH has been given in [104] by means of the computation of a topologic invariant, the Euler 
characteristic, in the case of a 2d lattice ip'^ model; (Hi) a direct analytic confirmation of the TH has been found in 
the particular case of the mean-field XY model [101]. We report on items (i), (ii) and (Hi) in §VIB 1, in §VIB2 and 
in §VIB3, respectively. 
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1. Indirect numerical investigations of the topology of configuration space 



In order to separate the singular effects due to the singular character of statistical measures at a phase transition 

from the singular effects due to some topological transition in configuration space, the first natural step is to consider 
again aj^ as observable, and to compute its integral value on suitable submanifolds of configuration space by means 
of a geometric measure, i.e. by means of a measure which has nothing to do with statistical ensemble measures. 

Consider as ambient space the iV-dimensional configuration space M of a Hamiltonian system with N degrees of 
freedom, which, when N oo, undergoes a phase transition at a certain finite temperature (or critical energy per 
degree of freedom e^), and let V — V{(p) be its potential energy per degree of freedom. 

Then the relevant geometrical objects are the submanifolds of M defined by 

Mu = V-\-oo,u] = {ip€M :V{ip)<u} , (123) 

i.e., each AIu is the set {(pi}fLi such that the potential energy per particle does not exceed a given value u. As u is 
increased from — oo to +00, this family covers successively the whole manifold M. All the submanifolds M„ can be 
given a Riemannian metric g whose choice is largely arbitrary. On all these manifolds {Mu,g) there is a standard 
invariant volume measure: 

dr] = ^/det{g) dq^ ...dq^ , (124) 

which has nothing to do with statistical measures. Let us finally define the hypersurfaces S„ as the u-level sets of V, 
i.e., 

^u = V-\u), (125) 

which arc nothing but the boundaries of the submanifolds M„. 

According to the discussion reported in § VI A, an indirect way to study the presence of topology changes in the 
family {(M„, (/)} is to look at the behaviour of the fluctuations of the Gaussian curvature, aj^, defined as 

= {Khh. - {KG)k (126) 

where (•) stands for integration over the surface as a function of u. The presence of cusp-like singularities of cr|- 
for some critical value of w, Uc, would eventually signal the presence of a topology change of the family {(M„, g)} at 
Uc- Such an indirect geometric probing of the presence of critical points seems an expedient way to probe the possible 
topology changes of the manifolds {Mu,g). In fact, the properties of the manifolds M„ are closely related to those of 
the hypersurfaces {Em}«<Uci can be inferred from the equation 

/ fdr,= r dv [ f\j:^d^/\\VV\\ (127) 

JM„ Jo JT,^ 

where dio is the induced measure"*^* on I]„ and / a generic function [105]. From Morse theory (see Appendix B) we 
know that the surface Yiu^ defined by V = Uc is a degenerate quadric, so that in its vic;inity some of the principal 
curvatures [103] of the surfaces Su~«c t^'^d to diverge^^. Such a divergence is generically detected by any function of 
the principal curvatures and thus, for practical computational reasons, instead of the Gauss curvature (which is the 
product of all the principal curvatures) we shall consider the total second variation of the scalar curvature TZ (i.e., the 
sum of all the possible products of two principal curvatures) of the manifolds {Mu,g), according to the definition 



cf\{u) = [Vol{Mu)]-^ ( dr] n-[Vol{Mu)]-' f drjTl 



(128) 



with TZ = g'^-' Rkij, where are the components of the Riemann curvature tensor [see Appendix A, Eq.(A39)] and 
Vol{Mu) = dr]. The subsets M„ of configuration space are given the structure of Riemannian manifolds {Mu,g) 



If a surface is parametrically defined tlirougli tlie equations = x'^{z^, . . . , z''), i = 1, . . . , N then the metric gij induced on 



the surface is given by g^j{z^, = YZ=i II^^ 



^^The principal curvatures are the inverse of the curvature radii measured, at any given point of a surface, in suitable directions. 
At a Morse critical point some of these curvature radii vanish. 
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by endowing all of them with the same metric tensor g. However, the choice of the metric g is arbitrary in view of 
probing possible effects of the topology on the geometry of these manifolds. 

In Ref. [99] the configuration spaces of a (p'^ model, defined on a one-dimensional lattice (linear chain) and on a 
two-dimensional square lattice have been considered. We recall that in the 2-d case this system undergoes a phase 
transition at a finite temperature, while in the 1-d case no phase transition is present. Three different types of metrics 
have been considered, i.e., 

(i) gj^J = [A — V{(p)]Sfi,y, i.e., a conformal deformation (see Appendix A, §A3) of the Euclidean flat metric 6^,^, 
where A > is an arbitrary constant chosen large enough to be sure that in the relevant interval of values of u 

the determinant of the metric is always positive; 

(2) (3) 

{ii) 9nv and giiv are generic metrics (no longer conformal deformations of the flat metric) defined by 

, A; = 2, 3 (129) 




where I is the N — 2 dimensional identity matrix, g^^^ is obtained by setting /'^^ = X^agz"* '■Pa + .9^'^' 
by setting Z'^' = ^ E aez<i fa + A with A> 0, and a labels the N lattice sites of a linear chain {d = 1) or of 
a square lattice {d = 2, N = n x n). 

These choices are completely arbitrary, however, and only if metrics of very simple form are chosen, both analytical 

and numerical computations are feasible also for rather large values of N. Thus the first metric has been chosen 
diagonal, and the other two metrics concentrate in only one matrix element all the non-trivial geometric information. 
Moreover, the first metric still contains a reference to the physical potential, whereas the other two define metric 
structures that are completely independent of the physical potential and only contain monomials of powers sufficiently 
high that they do not vanish after two successive derivatives have been taken (needed to compute curvatures). The 
topology of the subsets of points M„ and of is already determined (though well concealed) by the definitions of 
Eqs.(123) and (125); the task is to "capture" some information about their topology through a mathematical object 
or structure, defined on these sets of points, which is capable of mirroring the variations of topology through the 
u-pattern of an analytic function. This idea follows the philosophy of standard mathematical theories of differential 
topology, for example, within Morse theory, the information about topology is extracted through the critical points 
of any function - defined on a given manifold - fulfilling some conditions (necessary to be a good Morse function), 
or, whithin cohomology theory [106], topology is probed through vector spaces of differential forms (the de Rham's 
cohomology vector spaces) "attached" to a given manifold. Provided that good mathematical quantities are chosen 
as topology-variation detectors, arbitrary Riemannian metric structures could work as well. 

For the above defined metrics g'^''\ k = 1,2,3, simple algebra leads from the definition of the scalar curvature (see 
Appendix A) to the following explicit expressions: 



7e(i) = (iv - 1) 



Ay WVVP fN 3 



{A - vy (A - F)3 V 4 2 
||V/«|12 



(/(fc) - 1) 



2(/W - 1) 



(130) 

A; = 2, 3 (131) 



where V and A are the euclidean gradient and laplacian respectively, and V and A lack the derivative d/d^pa with 
a = 1 in the d—1 case, and lack the derivative d/dfa with a = (1, 1) in the d = 2 case. 

The numerical computation of the geometric integrals in Eq.(128) is worked out by means of a MonteCarlo algorithm 
[33,100] to sample the geometric measure drj by means of an "importance sampling" algorithm [107] suitably modified. 
In Figs. 23 and 24 (J^{u) are given for the one and two dimensional cases obtained for two different lattice sizes with 
^(1) (Fig. 23), and at given lattice size with 5(2,3) (pig_ 23). Peaks of cr^(u) appear at a certain value Vc, of Vc in the 
two-dimensional case, whereas only smooth patterns are found in the one- dimensional case, where no phase transition 
is present. 
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FIG. 23. Variance of the scalar curvature of Mu vs u computed with the metric ff'^'. Full circles correspond to the Id-cp* 
model with N = 400. Open circles refer to the 2d-ip^ model with TV = 20 x 20 lattice sites, and full triangles refer to 40 x 40 
lattice sites (whose values are rescaled for graphic reasons). 
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FIG. 24. <y^{u) of Mu vs u computed for the i^* model with: metric g^'^^ in Id, N = 400 (open triangles); metric g^^^ in 2d, 
AT = 20 X 20 (full triangles); metric 5^^^ in Id, N = 400 (open circles); metric g^^^ in 2d, TV = 20 x 20 (full circles). 



According to the discussion above, these peaks can be considered as indirect evidence of the presence of a topology 

transition in the manifolds M„ at u = Uc in the case of the two-dimensional (p^ model. It is in particular the 
persistence of cusp-like patterns of cr^{u) independently of the metric chosen that lends credit to the idea that this 
actually rcfic;c;ts a topological transition. Now we want to argue that the topological transition occurring at Uc is 
related to a thermodynamic phase transition which occurs in the model. In order to do that, in Ref. [99] the 
average potential energy per particle 

u{T) = (V) (132) 

has been numerically computed, as a function of T, by means of both MonteCarlo averaging with the canonical 
configurational measure, and Hamiltonian dynamics. In the latter case the temperature T is given by the average 
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kinetic energy per degree of freedom, and u is obtained as time average. Fig. 25 shows a perfect agreement between 
time and ensemble averages. The phase transition point is well visible at m = Mc — 3.75. Looking at Figs. 23 and 
24, we realize that, within the numerical accuracy, the critical value of the potential energy per particle Uc where 
the topological change occurs equals the statistical-mechanical average value of the potential energy at the phase 
transition. 




20 40 
T 

FIG. 25. Average potential energy vs temperature for the 2-d lattice tp"^ model with 0(1) symmetry. Lattice size AT = 20 x 20. 
The solid line is made out of 200 points obtained as time averages. Full circles represent MonteCarlo estimates of canonical 
ensemble averages. The dotted lines locate the phase transition. 

At this point the doubt, formulated at the beginning of this Section, about the possible non-geometrical origin 
of the "singular" , cusp- like patterns of {u) has been dissipated. These results have been found independently of 
statistical meclianical measures and of their singular cliarac;tc;r in presence of a phase transition. These results are also 
independent - at least to the limited extent of the three metric tensors reported above - of the geometric structure 
given to the family {M„}. Thus they seem most likely to have their origin at a deeper level than the geometric one, 
i.e. at the topologic level. Hence the observed phenomenology strongly hints that some major change in the topology 
of the configuration-space-submanifolds {M„} occurs when a second-order phase transition takes place. 

2. Direct numerical investigation of the topology of configuration space 

Though still based on numerical computations for a special model, a direct evidence of the tight relation between 
topology and phase transitions has been obtained by computing the M-dependence of a topologic invariant of the 
leaves E„ in the foliation of configuration space into a family of equipotential surfaces. 

In order to directly probe if and how the topology change - in the sense of a breaking of diffeomorphicity of the 
surfaces - is actually the counterpart of a phase transition, a diffeomorphism invariant has to be computed. This 
is a very challenging task because of the high dimensionality of the manifolds involved. Moreover, any algorithm of a 
combinatorial type (like those implied by simplicial decompositions, i.e. high dimensional analogs of tesselations with 
triangles that are used, for example, in numerical quantum gravity for low dimensional manifolds) is here hopeless. 
Only through a link between analytic and topologic mathematical objects can one hope to work out some direct 
information about topology. One such a link is provided by the Gauss-Bonnet-Hopf theorem that relates the Euler 
characteristic (see Appendix B) x(^«) with the total Gauss-Kronecker curvature of the manifold, i.e. [102] 

X(E„)=7^ Kg da (133) 

which is valid in general for even dimensional hypersurfaces of euclidean spaces [here dim(S„) = n = N — 1], 
and where: 7 = 2/yoZ(S") is twice the inverse of the volume of an n-dimensional sphere of unit radius; Kg is the 
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Gauss-Kroncckcr curvature of the manifold; da = ^ det{g)dx^dx^ ■ ■ ■ dx" is the invariant volume measure of E„ and 
g is the Riemannian metric induced from R''^. The Gauss-Kronecker curvature at a given point of a hypersurfacc is 
the product of the eigenvalues of its so-called shape operator; these eigenvalues are the principal curvatures of the 
hypersurface at the given point. The shape operator is constructed through the directional derivatives of the unit 
normal vector to the hypersurface at the given point computed in the n directions of the basis vectors of the plane 
tangent to the surface at the same point [103]. 

The numerical application of the Gauss-Bonnet-Hopf theorem has been worked out for a lattice (f'^ model in one 
and two spatial dimensions [100,112,104]. The main result is shown in Fig. 26 where xC^v) if^ reported vs v. In the Id 
case (open circles) a "smooth" pattern oi x{v) is found, whereas in the 2d case a cusp-like shaped x(u) shows up with 
the singular point corresponding to the phase transition point marked by the vertical dotted line. The parameters 
are those of the preceding Section, therefore the phase transition point is at Vc/N 3.75. 




v/N 

FIG. 26. Euler characteristic xi^v) for 1-d and 2-d tp^ lattice models. Open circles: 1-d case, N = 49; full circles: 2-d 
case, JV = 7 X 7. The vertical dotted line, computed separately for larger N , accurately locates the phase transition and the 
parameters are the same as in Figs. 23, 24 and 25. 

These results have two important consequences: i) the non constant value of x(w) in the Id case clearly shows that 
topology changes are there even in the absence of phase transitions; ii) an abrupt change in the rate of variation of 
topology with v seems the hallmark of a phase transition. Thus we have direct numerical evidence about the actual 
implication of topology in the appearance of a phase transition. At the same time we have evidence of a non simple 
one-to-one correspondence between topology changes and phase transitions. This is in full agreement with what is 
discussed in the next Section about a different kind of potential which can be analytically investigated. 



3. Topological origin of the phase transition in the mean-field XY model 

Until now we have not yet given any direct analytic evidence of the validity of the TH. Let us now consider again the 
mean- field XY model (116). In the case of this particular model it is possible to show analytically that a topological 
change in the configuration space exists and that it can be related to the thermodynamic phase transition [101]. 

Let us consider again, as was already done in §V1B1, the family of submanifolds of the configuration space 
defined in Eq. (123); now the potential energy per degree of freedom is that of the mean-field XY model, i.e., 

^("^^ = 4r = ^ ^ ~ ~ - /i^cos^i , (134) 

where G [0, 2tt\. Such a function can be considered a Morse function on M, so that, according to Morse theory (see 
Appendix B), all these manifolds have the same topology until a critical level V~^{vc) is crossed, where the topology 
of My changes. 
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A change in the topology of M^, can only occur when v passes through a critical value of V. Thus in order to detect 
topological changes in My we have to find the critical values of V, which means solving the equations 



d(Pi 



= 



N . 



(135) 



For a general potential energy function V, the solution of the Eqs. (135) would be a formidable task [108], but in the 
case of the mean- field XY model, the mean- field character of the interaction greatly simplifies the analysis, allowing 
an analytical treatment of the Eqs. (135); moreover, a projection of the configuration space onto a two-dimensional 
plane is possible. 

We recall that in the limit /i — > 0, the system has a continuous phase transition, with classical critical exponents, 
at Tc = 1/2, or Sc = 3/4, where e = E/N is the energy per particle. We aim at showing that this phase transition 
has its foundation in a basic topological change that occurs in the configuration space M of the system. Let us 
remark that since V{(p) is bounded, —h < V(lp) < 1/2 + h'^/2, the manifold is empty as long as v < —h, and when v 
increases beyond 1/2 -t- no changes in its topology can occur so that the manifold My remains the same for any 
V > 1/2 -|- h^/2, and is then an A^-torus. To detect topological changes we have to solve Eqs. (135). To this end it is 
useful to define the magnetization vector, i.e., the collective spin vector m = Si> which as a function of the 

angles is given by 



Due to the mean-field character of the model, the potential energy (116) can be expressed as a function of m alone 
(remember that J — 1), so that the potential energy per particle reads 



Now we can solve these equations and find all the critical values of V. The solutions of Eqs. (138) can be grouped in 
three classes: 

(i) The minimal energy configuration ift = Vz, with a critical value v = Vq = —h, which tends to as /i ^ 0. In 
this case, -I- = 1. 

(ii) Configurations such that my — 0,smipi — Vi. These are the configurations in which ipi equals either or 
tt; i.e., we have again ip^ = Vi, but also the configurations with (p^. = n and (fi = k, as well as the 
N{N—1) configurations with 2 angles equal to tt and all the others equal to 0, and so on, up to the configuration with 

= TT Vi. The critical values corresponding to these critical points depend only on the number of tt's, n^, so that 
v{nTr) = ^[1 — ^{N — 2n^)^] — ;^(A^ — 271^^)- We see that the largest critical value is, for A^ even, v(n7r = A^/2) = 1/2 
and that the number of critical points corresponding to it is 0{2^). 

(Hi) Configurations such that nix = —h and niy = 0, which correspond to the critical value Vc = 1/2 + h'^/2, which 
tends to 1/2 as /i ^ 0. The number of these configurations grows with A^ not slower than A^! [101]. 

Configurations (i) are the absolute minima of V, (iii) are the absolute maxima, and (ii) are all the other stationary 
configurations of V. 

Since for v < vq the manifold is empty, the topological change that occurs at vq is the one corresponding to the 
"birth" of the manifold from the empty set; subsequently there are many topological changes at values w(n^) € (vq, 1 /2] 
till at Vc there is a final topological change which corresponds to the "completion" of the manifold. We remark that 
the number of critical values in the interval [vq, 1/2] grows with N and that eventually the set of these critical values 
becomes dense in the limit A^ oo. However, the critical value Vc remains isolated from other critical values also in 
that limit. We observe that it is necessary to consider a nonzero external field h in order that V is a Morse function, 
because if h = all the critical points of classes (i) and {ii) are degenerate, in which case topological changes do 
not necessarily occur^°. This degeneracy is due to the 0(2)-invariance of the potential energy in the absence of an 



It would also be possible to avoid this problem by considering an improved version of Morse theory, referred to as equivariant 
Morse theory [109]. 




(136) 



V{ip) =V{m^,my) = -{1 - ml - ml) - hm, 



(137) 



This allows us to write the Eqs. (135) in the form {i = 1, . . . ,N) 

{rrix + h) sin ipi — niy cos </3j = . 



(138) 
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external field. To be sure, for /i ^ 0, V may not be a Morse function on the whole of M cither, but only on with 
V < Vc, because the critical points of class {Hi) may also be degenerate, so that Vc does not necessarily correspond to 
a topological change. However, this difficulty could be dealt with by using that the potential energy can be written in 
terms of the collective variables and my — as in Eq. (137). This implies that we consider the system of N spins 
projected onto the two-dimensional configuration space of the collective spin variables. According to the definition 
(136) of m, the accessible configuration space is now not the whole plane, but only the disk 

D = {{m:„my):ml+ml<l} . (139) 

Thus we want to study the topology of the submanifolds 

Dy = {{mx,my) & D ■.V{mx,my) <v} . (140) 




V = < i; < 1/2 17 = 1/2 



FIG. 27. The sequence of topological changes undergone by the manifolds Dv with increasing v in the limit /i — > 0. 

The sequence of topological transformations undergone by Dy can now be very simply determined in the limit /i ^ 
(see Fig. 27), as follows. As long as w < 0, Dy is the empty set. The first topological change occurs at w = uq = 
where the manifold appears as the circle + = 1, i.e., the boundary dD of D. Then as v grows Dy is given by 
the conditions 

l-2v <ml+ ml <1 , (141) 

i.e., it is the ring with a hole centered at (0,0) (punctuated disk) comprised between two circles of radii 1 and 
respectively. As v continues to grow the hole shrinks and is eventually completely filled as u = Vc = 1/2, where the 
second topological change occurs. In this coarse-grained two-dimensional description in D, all the topological changes 
that occur in M between v = and v = 1/2 disappear, and only the two topological changes corresponding to the 
extrema of V, occurring a,t v = vq and v = Vc, survive. This means that the topological change at Vc should be present 
also in the full A'^-dimensional configuration space, so that the degeneracies mentioned above for the critical points of 
class {Hi) should not prevent a topological change. 

Now we want to argue that the topological change occurring at Vc is related to the thermodynamic phase transition 
of the mean-field XY model. Since the Hamiltonian is of the standard form (1), the temperature T, the energy per 
particle e and the average potential energy per particle u = (V) obey, in the thermodynamic limit, the following 
equation: 

e=^+u{T), (142) 

where we have set Boltzmann's constant equal to 1. Substituting the values of the critical energy per particle Ec = 3/4 
and of the critical temperature = 1/2 we get = u{Tc) — 1/2, so that the critical value of the potential energy per 
particle Vc where the last topological change occurs equals the statistical-mechanical average value of the potential 
energy at the phase transition, 

Vc = Uc . (143) 

Thus although a topological change in M occurs at any N, and Vc is independent of A'^, there is a connection of 
such a topological change and a thermodynamic phase transition only in the limit A ^ oo, ft ^ 0+, when indeed a 
thermodynamic phase transition can be defined. 

A similar kind of difference, as here between topological changes in mathematics (for all N) and phase transitions 
in physics (for N oo only), also occurs in other contexts in statistical mechanics, e.g. in nonequilibrium stationary 
states [110]. 
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Since not all topological changes correspond to phase transitions, those that do correspond, remain to be determined 
to make the conjecture of Ref. [36] more precise. In this context, we consider one example where there are topological 
changes very similar to the ones of our model but no phase transitions, i.e., the one-dimensional XY model with 
nearest-neighbor interactions, whose Hamiltonian is of the class (1) with interaction potential 

^ N N 

^(^) = 4 23 " ^08(^1 -^i)]-hY^ COS . (144) 

i=l i=l 

In this case the configuration space M is still an A/'-torus, and using again the potential energy per degree of freedom 

V = V/N as a Morse frmction, we can see that also here there are many topological changes in the submanifolds 
My as V is varied in the interval [0,1/2] (after taking h 0+). However there arc critical points of the type 
ipj = ipk = ipi = ■ ■ ■ = TT, ipi = yi ^ i, k,l, . . .; at variance with the mean-field XY model, it is now no longer 
the number of tt's that determines the value of V at the critical point, but rather the number of domain walls, n^, 
i.e., the number of boundaries between "islands" of tt's and "islands" of O's: v{nd) = Ud/^N. Since G [0, A''], the 
critical values lie in the same interval as in the case of the mean-field XY model; but now the maximum critical value 

V = 1/2, instead of corresponding to a huge number of critical points, which rapidly grows with N, corresponds to 
only two configurations with N domain walls, which are </?2fe = 0, ip2k+i = tt, with A; = 1, . . . , N/2, and the reversed 
one. There are also "spin-wavc-likc" critical points, i.e., such that 6*^*= = const e27rifcn/Ar ^ — l,...,iV [111]; 
their critical energies are comprised in the interval above but again there is not a critical value associated to a huge 
number of critical points. 

Thus this example suggests the conjecture that a topological change in the configuration space submanifolds Aly 
occurring at a critical value Wc, is associated with a phase transition in the thermodynamic limit only if the number 
of critical points corresponding to the critical value Vc is sufficiently rapidly growing with N . On the basis of the 
behavior of the mean-field XY model we expect then that such a growth should be at least exponential. Furthermore, 
a relevant feature appears to be that Vc remains an isolated critical value also in the limit N ^ oo: in the mean- field 
XY model this holds only if the thermodynamic limit is taken before the /i 0+ limit: this appears as a topological 
counterpart of the non-commutativity of the limits /i ^ 0+ and iV ^ oo in order to get a phase transition in statistical 
mechanics. 

The sequence of topological changes occurring with growing V makes the configuration space larger and larger, till 
at Vc the whole configuration space becomes fully accessible to the system through the last topological change. From 
a physical point of view, this corresponds to the appearance of more and more disordered configurations as T grows, 
which ultimately lead to the phase transition at Tc. We remark that the connection between the topology of the 
configuration space and the physics of continuous phase transitions made here via the potential energy, in particular 
Eq. (143), only makes sense in the thermodynamic limit, where the potential energy per particle u{T) is well-defined 
since its fluctuations vanish then at least as l/v^- This holds for our mean-field model, since for such a model, all 
fiuctuations are absent. In the case of a real continuous (critical) phase transition the non-trivial role of fluctuations 
may complicate the present picture. 



C. The Topological Hypothesis 

The statistical behaviour of physical systems described by standard Hamiltonians is encompassed, in the canonical 
ensemble, by the partition function in phase space 




(145) 



where p= (pi, . . . ,pjv), dp = lli=idPi, Q = (gi, • . • , q-iv), dq = Ui^idqi; = {(gi, . • . , Q-jv) e R^|y(gi, . . . , gjv) = u} 
are equipotcntial hypcrsurfaccs of configuration space and da is the measure on S„ induced from R^. The last 
term of the equation above shows that for Hamiltonians (1) the relevant statistical information is contained in the 
canonical configurational partition function = J dqexp[—l3V{q)]. Moreover, the last term of Eq. (145), written 
using a co-area formula [105], decomposes Z^jj into an infinite summation of geometric integrals da /||Vy||. This 
decomposition provides the point of attack for the formulation of a general hypothesis about a deep link between 
geometry, topology and thermodynamics and, obviously, phase transitions. In fact, once the potential energy V{q) is 
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given, the configuration space of the system is antomatically fohated into the family {Euj^gR of these cquipotcntial 
hypersurfaces. Now, from standard statistical mechanical arguments concerning the equivalence of canonical and 
microcanonical ensembles we know that, at any given value of the inverse temperature /3, the larger the number 
N of particles, the closer to some are the microstates that significantly contribute to the statistical averages of 
thermodynamic observables. The hypersurface is the one associated with u = Uf} = (Z^)^^ J JJ dqiV{q)e~l^^^'^\ 
the average potential energy computed at a given /3. Thus, at any /?, if N is very large the effective support of the 
canonical measure shrinks very close to a single = E„^. Hence and on the basis of what was found in [36,101,99], 
it was formulated the following 

Topological Hypothesis: The basic origin of a phase transition lies in a topological change of the support of the measure 
describing a system. This change of topology induces a, change of the m,easure itself at the transition point. 

In other words, this hypothesis stipulates that some change of the topology of the {E„}, occurring at some = 
Uc{Pc), could be the origin of the singular behavior of thermodynamic observables at a phase transition rather than 
measure singularities which in this view are induced from a deeper level where the topology changes take place. 

In other words, the claim is that the canonical measure should "feel" a big and sudden change - if any - of the 
topology of the equipotential hypersurfaces of its underlying support, with as a consequence, the appearance of the 
typical signals of a phase transition, i.e. almost singular (at finite N) energy or temperature dependences of the 
averages of appropriate observables. The larger TV, the narrower the effective support is of the measure - as discussed 
above - and hence the sharper the mentioned signals can be, until true singularities appear in the N oo limit. 

We emphasize though that not all topological transitions lead to physical phase transitions. At present the precise 
connection between topological transitions and phase transitions still has to be clarified in many respects. Certain is 
that not every topological transition corresponds to a phase transition, as has been discussed in Sections VI B 2 and 
VI B 3. Rather it seems that, on the basis of present evidence, a phase transition corresponds to a super-combination of 
many simultaneous topological transitions taking place, where many might mean at least exponentially growing with 
the number of particles. It seems therefore more like a super-topologically constructed transition. This is illustrated 
analytically by the above discussed XY mean-field model, where an exponential crowding of topological transitions 
occurs on one side of the phase transition. Though such an analysis has not been possible for a numerically treated 
lattice if^ model, on the other hand, like the Euler characteristic x clearly shows, a phase transition corresponds to 
an abrupt transition between different rates of change in the topology above and below of the phase transition point; 
no such information is available for the analj^tic XY mean-field model where the calculation of x is for the moment 
very difficult. 

D. Open questions and future developments 
1. The phase space trajectories of dynamical systems described by Hamiltonian functions of the form (1), i.e., 

1 ^ 

^ i=l 

where the t^'s and the tt's are, respectively, the coordinates and the conjugate momenta, are bound to the constant 
energy hypersurfaces Ee of the 2iV-dimensional phase space spanned by the i^'s and the tt's. Therefore it would 
be natural to investigate the relationship between the topology variations of these hypersurfaces in phase space and 
phase transitions. The Hamiltonian H would be the Morse function in this case. Moreover, having already considered 
the role played by the S„ of configuration space, we can wonder what is the relationship between topology changes in 
configuration space and phase space respectively. Such a relationship is somehow subtle. At first sight all the critical 
points of configuration space are embodied also in critical points of phase space, in fact all the critical points of H are 
such that, Vi, tt^ = and \7,-V{ip) = 0. However, the critical points of 7Y physically correspond to vanishing kinetic 
energy, so that, if those topology changes of the S„ that are associated to a phase transition were to correspond to 
topology changes of the S^; - because the critical points of H incorporate those of V{(p) - then the critical potential 
energy density and the critical total energy density at a phase transition should coincide, which is not the case. 
An argument to clarify this point can be given [50] by observing that the dynamics does not equally sample the 
whole Eg: the larger N, the smaller the relative fluctuations of the potential energy {5^Vy^^/{V) and kinetic energy 
{5^K)^/^ /{K), respectively, are. Thus, by putting v = {V) and t = {K), we can assume that, at large and given 
energy E, the momenta mainly live close to the hypersphere S^~^ = {(pi, . . . ,Pn) G X^iXi \Pi — ^"^^ 
lagrangian coordinates mainly live close to the equipotential hypersurface E^~^ with v + t = E. Therefore, though 
the microcanonical measure mathematically extends over a whole energy surface, as far as physics is concerned, at 
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very large A'^ a non-negligible contribution to the microcanonical measure is in practice given only by a small subset 
of an energy surface. This subset can be approximately modeled by the product manifold Sf-i X "\ Since the 
kinetic energy submanifold Sf ~^ is a hypersphere at any t, a change in the topology of S^~^ directly entails a change 
of the topology of x Sf'"^, that is of the effective model-manifold for the subset of where the dynamics 

mainly "lives" at a given energy E. 

This confirms that, as long as we are interested only in classical Hamiltonian systems of the standard form (1), 
we can restrict our geometrical and topological investigation to the submanifolds of the A''-dimensional configuration 
space M. 

2. A recent and important advance has been achieved at a more mathematical level. In Refs. [100,112,113] a 
theorem has been proved that establishes the necessity of topological changes of the equipotential hypersurfaces S„ 
for the appearance of first or second order thermodynamic phase transitions. The theorem applies to a wide class of 
finite-range potentials, bound below, describing systems confined in finite regions of space with continuously varying 
coordinates. The proof proceeds by showing that, under the crucial assumption of diffeomorphicity of the I]„ in an 
arbitrary interval of values for u, the Helmoltz free energy is uniformly convergent in N to its thermodynamic limit, 
at least within the class of twice differentiablc functions, in a corresponding interval of temperature. 

This theorem confirms the general validity of the TH and ensures that for a wide class of physical potentials 
the mathematical framework of differential topology is adequate to describe, at least, first and second order phase 
transitions. There is no proof of sufficiency. On the basis of the discussions in Sections VI B 2, VI B 3 and VIC we 
already know that a simple loss of diffeomorphicity of the S„ is not sufficient to lead to a phase transition. 

A thorough investigation of those classes of topology changes that are responsible for the appearance of phase 
transitions is at present the main challenge of this new point of view about phase transitions and certainly represents 
a topic that will remunerate the efforts addressed to it. 

3. Let us finally highlight some interesting related topics. 

(i) We might speculate about the possibility of relating universal quantities of the theory of critical phenomena, 
like critical exponents, to some topological counterpart; in fact a notion of universality arises quite naturally in 

a topological framework. 

(a) Topology provides a common ground for the roots of both dynamics and thermodynamics: insofar as the 
dynamics of the system, i.e. the motion of the trajectory in phase space, takes place in what was called before the 
support of the statistical measure, it is clear that the nature of the trajectory will crucially depend on the topology 
of the manifold to which it belongs. This therefore strenghtens the interest of a dynamical treatment of phase 
transitions giving new emphasis to the microcanonical ensemble and thereby joining other recent developments in 
the field [114]. Since the dynamical approach does not depend on whether a system is in statistical equilibrium, 
non-equilibrium, or in a metastable state (like a glass and, more generally, amorphous materials) a dynamical 
approach to phase transitions might also be important for systems whose thermodynamical state is not well 
defined. 

(iii) The fact that the topological changes appear at any N opens a new possibility to study transitional phenomena 
in finite systems, like nuclear and atomic clusters, polymers and proteins, or other biological systems, as well as 
for nano and mesoscopic structures. 

ACKNOWLEDGMENTS 

We want to dedicate this paper to the memory of our friend and brilliant collaborator Lando Caiani, prematurely 
deceased, who gave seminal contributions to many aspects of the approach described here, and especially to the early 
development of the connection between topology and phase transitions. 

During their development, the main ideas described in the present Report have been discussed with many colleagues, 
among which we would like to thank, for their particular interest and for many useful suggestions and criticism, and 
sometimes for an active collaboration, S. Caracciolo, M. Ccrruti-Sola, C. Clemcnti, Y. Elskcns, M.-C. Firpo, R. 
Franzosi, R. Gatto, R. Livi, G. Mussardo, H. A. Posch, V. Penna, G. Pettini, M. Rasetti, S. Ruffo, L. Spinelli, G. 
Vezzosi, A. M. Vinogradov, and the DOCS research group in Firenze (http://docs.de. unifi.it/^docs/). 

Part of this work has been done at the Erwin Schrodinger Institut fiir Mathematische Physik (ESI) in Vienna, 
whose kind hospitality is gratefully acknowledged. LC thanks the Rockefeller University in New York for its kind 
hospitality and for partial financial support, the 1997-1999 SINTESI program of the Italian Ministry of University and 
of Scientific and Technologic Research (MURST) for partial financial support, and Prof. H. A. Posch for hospitality 
at the Institut fiir Experimental Physik of the University of Vienna. EGDC is indebted to the Basic Engineering 



55 



Section of the US Department of Energy for support under grant DE-FG02-88-ER13847, as well as to the Politecnico 
di Torino and to the INFM, UdR Torino Politecnico, for hospitality and partial financial support. 



APPENDIX A: SUMMARY OF RIEMANNIAN GEOMETRY 

In the following we briefly recall some essential concepts and notations of Riemannian differential geometry which 
are used in the main text. The present section is only meant to facilitate the reader to follow the main text of the 
Report, so that our discussion will not be a rigorous treatment of the subject. For a more elaborate discussion, we 
refer the reader to a textbook of general relativity (e.g., the classic Landau and Lifshitz's book [115] or the more 
recent and complete, but still very clear and readable, textbook by Wald [116]). A more mathematically oriented, 
but still readable by physicists, introduction to the subject is given by do Carmo [58]; a comprehensive and rigorous 
treatment, which, however, goes far beyond what is needed to follow the exposition in the main text, can be found in 
Kobayashi and Nomizu [117]. 

The Einstein summation convention over repeated indices is always understood unless explicitly stated to the con- 
trary. Moreover, we follow throughout the paper the usual convention to suppress the dependence of the components 
of vector and tensor quantities on the (proper) time and, in general, only indicate it explicitly when this dependence 
is absolutely relevant. 

1. Riemannian manifolds 

A set M is called a differentiable manifold if it can be covered with a collection, cither finite or denumerable, of 
charts, such that each point of M is represented at least on one chart, and the different charts are differentiably 
connected to each other. A chart is a set of coordinates on the manifold, i.e., it is a set of n real numbers {x^,. . . , x") 
which denote the "position" of a point on the manifold. The number n of coordinates of a chart is the same for each 
connected part of the manifold (and for the whole manifold if the latter is connected, i.e., it cannot be split in two 
disjoint parts which are still manifolds); such a number is called the dimension of the manifold M. The union of the 
charts on M is called an atlas of M. 

a. Vectors and tensors 

A vector (more precisely, a tangent vector), can be defined using curves on the manifold M. Given a curve 7 in M, 
represented in local coordinates by the parametric equations x = ip{t), we define a tangent vector at P e M as the 
velocity vector of the curve in P, i.e., 

. = 7 = lim^^(^l^, m = P, (Al) 
so that the n components of the tangent vector v are given by 

The set of all the tangent vectors of M in P is a linear space, referred to as the tangent space of M in P, and denoted 
by TpM. Each tangent space is isomorphic to an n-dimensional Euclidean space. Given a chart (x^, . . . ,a;") in a 
neighborhood of P, a basis {Xi, . . . , X„) of TpM can be defined, so that a generic vector v is expressed as a sum of 
the Xi's weighted by its components, 

V = v'Xi . (A3) 

The basis {Xi} is called a coordinate basis of TpM, and its components Xi are often denoted^^ by d/dx^. The basis 
depends on the chart: choosing another chart, {x'^ , . . . ,x'"'), we get another basis {X'i}. The components of v in the 
two different bases are connected by the following rule, 



'^The origin of this notation is in the fact that vectors can be defined as directional derivatives on M (see e.g. Ref. [116]). 
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referred to as the vector transformation rule. Indeed, one can define a vector as a quantity whose components 
transform according to Eq. (A4). The union of all the tangent spaces of the manifold M, 



TM = y TpM , (A5) 

is a 2n-dimensional manifold and is referred to as the tangent bundle of M. 

A vector field 1/ on M is an assignment of a vector at each point P e M. If / is a smooth function, 

V{f)\p = vp{f) (A6) 

is a rcial number for each P S M, i.e., v{f) is a function on M. If such a function is smooth, V is called a smooth 
vector field on M. The curves ip{t) which satify the differential equations 

^ = V{^{t)) (A7) 

arc called the trajectories of the field V, and the mapping ipt : M M which maps any point P of M along the 
trajectory of V emanating from P is called the flow of V. Given two vector fields V, W, one can define the commutator 
as the vector field [V, W] such that 

[V, W]if) = V{W{f)) - WiVif)) , (A8) 

i.e., in terms of the local components, 

We note that, if {Xi} is a coordinate basis, 

[Xi,X,]=0 Vi,j, (AlO) 

and that, conversely, given n nonvanishing and commuting vector fields which are linearly independent, there always 
exists a chart for which these vector fields are a coordinate basis. 

Tangent vectors are not the only vector-like quantities that can be defined on a manifold M: there arc also cotangent 
vectors, which can be defined as follows. Let us recall that the dual space V* of a vector space V is the space of linear 
maps from V to the real numbers. Given a basis of V, {w,}, a basis of V*, {w**}, called the dual basis, is defined by 

u"{uj) =5] . (All) 

The dual space of TM, T*M, is called the cotangent bundle of M. Its elements are called cotangent vectors, or 
sometimes covariant vectors (while the tangent vectors are sometimes denoted as contravariant vectors). The dual 
basis elements are usually denoted as dx^ , . . . ,dx'", i.e., dx^ is such that dx^{d/dx^) = 6j. The components of 
cotangent vectors transform according to the rule 



dx^ 



to be compared with Eq. (A4). The common rule is to use subscripts to denote the components of dual vectors and 

superscripts for those of vectors. 

A {k, lytensor T over a vector space F is a multilinear map 

T -.V* ■■ ■ xV* xV X ■■■xV , (A13) 

^ V ' ^ V ' 

k times / times 

i.e., acting on k dual vectors and / vectors, T yields a number, and it does so in such a manner that if we fix all but 
one of the vectors or dual vectors, it is a linear map in the remaining variable. A (0,0) tensor is a scalar, a (0, 1) 
tensor is a vector, and a (1, 0) tensor is a dual vector. The space T{k, I) of the tensors of type (fc, I) is a linear space; 
a (fc, i)-tensor is defined once its action on k vectors of the dual basis and on I vectors of the basis is known, and 
since there are n'^n' independent ways of choosing these basis vectors, T{k, I) is a n'^^' -dimensional linear space. Two 
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natural operations can be defined on tensors. The first one is called contraction with respect to the i-th (dual vector) 
and the j-th (vector) arguments and is a map 

C :T GT{k,l)^CT &T{k-l,l-l) (A14) 

defined by 

n 

CT = ^T(...,^,...;...,^,...) . (A15) 

The contracted tensor CT is independent of the choice of the basis, so that the contraction is a well-defined, invariant, 

operation. The second operation is the tensor product, which maps an element T{k, I) x T{k' , I') into an element of 
T(fc + k' ,1 + 1'), i.e., two tensors T and T' into a new tensor, denoted by T ig) T' , defined as follows: given k + k' dual 
vectors v^*, . . . , * and I + I' vectors wi, . . . , , then 

T ® T'{v'* v''+'''*; wu..., wi+r ) = T{v^\ . . . y*;wu ■ ■ ■ ,wi) T'{v''+'\ . . . , v''+^'*; wi+,, wi+r ) . (A16) 

The tensor product allows one to construct a basis for T{k, I) starting from a basis {v^} of V and its dual basis {v'^*}: 
such a basis is given by the n'^"'"' tensors {Vfj,^ ® ■ ■ - ^ v^i^ ® v'^^* ® ■ ■ ■ ® v'^'-*} . Thus, every tensor T e T(fc, I) allows a 
decomposition 

n 

/ii,...,i'i=i 

the numbers T'*^'"'"'^^...!,, are called the components of T in the basis {v^}. The components of the contracted tensor 
CT are 

(remember the summation convention) , and the components of the tensor product T ®T' are 

(r®r)'^"''^'=+'=;,....,+,, =^'^-•''^,....,^''"=+^•••'''=+'=^,^,....,^,, . (ai9) 

All these results are valid for a generic vector space, so that they hold in particular for the vector spaces of the tangent 
bundle TM of M, over which tensors (and, analogously to vector fields, tensor fields) can be defined exactly as above. 

h. Riemannian metrics 

The infinitesimal square distance on Ad , i.e., the length element ds"^ (also referred to as the metric)) can be defined 
at each point P S M by means of a (0, 2)-tensor provided it is symmetric, i.e., g{v, w) = g{w, v), and nondegenerate, 
i.e., g{v,w) = Vw e TpM if and only ii w = 0^^. In fact, a g with these properties induces on the tangent bundle 
TM a nondegenerate quadratic form (called the scalar product), 

g{v, w) = {v, w) : TM x TM ^ R . (A20) 

Then it is possible to measure lengths on the manifold. A manifold M, equipped with a scalar product, is called a 
(psendo) Riemannian manifold, and the scalar product is referred to as a (paciido) Riem.annian structure on M . If the 
quadratic form (A20) is positive-definite, then one speaks of a (proper) Riemannian metric. In the latter case the 
squared length element is always positive. For instance, one can define the length of a curve as 



^(7) = / V(Mdt . (A21) 



^^Or, equivalently, g{v,w) = Vw e TpM if and only if v = 0; the two statements are equivalent because g is a symmetric 
tensor. 
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The curves 7 which are extremals of the length functional are called the geodesies of M. 
In a coordinate basis, we can expand the metric g as 

g = gij dx" ® dx' , (A22) 
so that one defines the invariant (squared) length element on the manifold, in local coordinates, as 

rfs2 = g.jdx'dx^ . (A23) 
The scalar product of two vectors v and w is given, in terms of g, by 

{VjW) = gijV^w^ = VjW^ = vhvi . (A24) 

In the above equation we have made use of the fact that g estabilishes a one-to-one correspondence between vectors 
and dual vectors, i.e., in components, 

gijV^ = Vi . (A25) 

For this reason, the components of the inverse metric g~^ are simply denoted by g'^\ instead of (fif"^)*-', and allow to 
pass from dual vector (covariant) components to vector (contravariant) components: 

g'^Vj = v' . (A26) 

This operation of raising and lowering the indices can be applied not only to vector, but also to tensor components. 
This allows us to pass from {k, I) tensor components to the corrseponding {k + 1,1 — 1) tensor components and vice 
versa. What does not change in the operation is the sum k + 1 which is called the rank (or the order) of the tensor. 



2. Covariant differentiation 



The introduction of a differential calculus on a manifold which is not Euclidean is complicated by the fact that 

ordinary derivatives do not map vectors into vectors, i.e., the ordinary derivatives of the components of a vector w, 
dw^ /dt, taken for instance at a point P along a given curve 7(t), are not the components of a vector in TpM , because 
they do not transform according to the rule (A4). The geometric origin of this fact is that the parallel transport of a 
vector from a point P to a point Q on a non-Euclidean manifold does depends on the path chosen to join P and Q. 
Since in order to define the derivative of a vector at P, we have to move that vector from P to a neighboring point 
along a curve and then parallel-transport it back to the original point in order to measure the difference, we need a 
definition of parallel transport to define a derivative; conversely, given a (consistent) derivative, i.e., a derivative which 
maps vectors into vectors, one could define the parallel transport by imposing that a vector is parallel transported 
along a curve if its derivative along the curve is zero. The two ways are conceptually equivalent: we follow the first 
way, by introducing the notion of a connection and then using it to define the derivative operator. Such a derivative 
will be referred to as the covariant derivative. 

A (linear) connection on the manifold M is a map V such that, given two vector fields^^ A and B, it yields a third 
field V aB with the following properties: 

1. WaB is bilinear in A and B, i.e., V A{aB + I3C) = aV aB + 13V aC and V aA+fisC = aV aC + fiV bC; 

2. Vy(^)B = /(V^B); 

3. (Leibnitz rule) W Af{B) = {dAf)B + /{VaB), where Oa is the ordinary directional derivative in the direction 
of A. 

The parallel transport of a vector V along a curve 7, whose tangent vector field is 7, is then defined as the (unique) 
vector field W{t) = W{'y{t)) along 7(f) such that 

1. W{0) = V; 



One could also consider tensor fields, but for the sake of simplicity we define connections using vectors. 
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V^V . (A27) 



2. V-yW = along 7. 

The notion of covariant derivative now immediately follows: the covariant derivative DV/ dt of V along 7 is given by 
the vector field 

DV 

~dr 

On the basis of Eq. (A27), with a certain abuse of language, one often refers to Vx^ as the covariant derivative of 
Y along X, where X and Y are generic vector fields. Among all the possible linear connections, and given a metric 
g, there is one and only one which (i) is symmetric, i.e., 

V^y - VyX = [X, Y] \/X, Y , (A28) 

and (a) conserves the scalar product, i.e., the scalar product of two parallel vector fields P and P' is constant along 

7, 

|(P,P') = 0. (A29) 

Such a linear connection is obviously the natural one on a Riemannian manifold, and is referred to as the Levi-Civita 
(or Riemannian) connection. Whenever we refer to a covariant derivative without any specification, we mean the 
covariant derivative induced by the Riemannian connection. 

The components of the Riemannian connection V with respect to a coordinate basis {Xi} are the Christoffel 
symbols, given by 

T% = {dq\Vx,Xk) . (A30) 
and are given, in terms of the derivatives of the components of the metric, by the following formula 

Tj-fe = ^ff"" [djgkm + dkgmj - dmgjk) , (A31) 

where di = d/dx^. The expression in local coordinates of the covariant derivative (A27) of a vector field V is then 

DV^ dV^ ■ dr^ , 



a. Geodesies 

The geodesies, which were already defined as the curves of extremal length on the manifold, can also be defined as 
self-parallel curves, i.e., curves such that the tangent vector 7 is always parallel transported. Thus geodesies are the 
curves 7(f) which satisiy the equation (referred to as the geodesic equation) 

f=0 (A33) 

whose expression in local coordinates foUws from Eq. (A32), and is 

rf^x' . dx^ dx'' „ 

Since the norm of the tangent vector 7 of a geodesic is constant, Idj/dt] = c, the arc lenght of a geodesic is proportional 
to the parameter: 



sit) 



rt2 
Jti 



d'j 
'dt 



dt = c{t2 -tl) . (A35) 



When the parameter is actually the arc length, i.e., c = 1, we say that the geodesic is normalized. Whenever we 
consider a geodesic, we assume it is normalized, if not explicitly stated otherwise. This means that the Eqs. (A34) 
are nothing but the Euler-Lagrange equations for the length functional along a curve 7(5) parametrized by the arc 
length, 

^(7) = [ ds . (A36) 

Given a geodesic 7(5) on M, there exists a unique vector field G on TM such that its trajectories are (7(5), 7(5)). 
Such a vector field is called the geodesic field and its flow the geodesic flow on M. 
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3. Curvature 



The curvature of a Riemannian manifold (M, g) is — intuitively — a way of measuring how much the manifold 
deviates from being Euclidean. The curvature tensor, also known as the Riemann-Christoffel tensor, is a tensor of 
order 4 defined as 

R{X, Y) = Vx Vy - VyVx - y[x,Y] , (A37) 

where V is the Levi-Civita connection of M. Observe that if M = R^, then R{X,Y) = for all the pairs of 
tangent vectors X, Y, because of the commutativity of the ordinary derivatives. In addition, R measures the non- 
commutativity of the covariant derivative: in fact, if we choose a coordinate system {xi, . . . ,Xn}, we have, since 
-9- -2-1=0 

dxi ' dxj ' 

f d d \ 

^[dxl'dx]) ^ ^9/3^* ~ Va/a^^.Va/axi • (A38) 
In local coordinates, the components of the Riemann curvature tensor (considered here as a (l,3)-tensor) are given by 

^'jkl = ~ + '^'jl^lr - '^kl^)r ■ (A39) 

Thus, given a metric g, the curvature R is uniquely defined. A manifold (M, g) is called flat when the curvature tensor 
vanishes. 

Given a positive function a conformal transformation is the transformation 

{M,g)^{M,g); g = fg . (A40) 

Two Riemannian manifolds are said conformally related if they are linked by a conformal transformation. In particular, 

a manifold is (M, g) conformally fl,at if it is possible^ to find a conformal transformation which sends g into a flat 
metric. Conformally flat manifolds exhibit sonic remarkable simplifications for the calculation of the curvature tensor 
components (sec e.g. [63]; an application is given in §111 C). 

Closely related to the curvature tensor is the sectional — or Riemannian — curvature, which we define now. Let 
us consider two vectors u,v G TpM, and let us put 

\uAv\ = {\u\'\v\^-{u,v)Y^' , (A41) 

which is the area of the two-dimensional parallelogram determined by u and v. If |w, A v\ ^ the vectors u, v span a 
two-dimensional subspace tt C TpM. We define the sectional curvature at the point P relative to tt, as the quantity: 

i.(P;.,.) = K(P,.)=<5Ml^ 

which can be shown to be independent of the choice of the two vectors u,v G tt. In local coordinates, Eq. (A42) 
becomes 

K{P; u, v) = Rijki , , ,2 • (A43) 

The knowledge of K for the A''(A'' — 1) planes tt spanned by a maximal set of linearly independent vectors completely 
determines R at P. 

If dim(M) — 2 then K coincides with the Gaussian curvature of the surface, i.e., with the product of the reciprocals 
of two curvature radii. 

A manifold is called isotropic if K{P, tt) does not depend on the choice of the plane tt. The remarkable result — 
Schur's theorem [58] — is that in this case K is also constant, i.e. it does not depend on the point P either. 

Some "averages" of the sectional curvatures arc very important. The Ricci curvature Kp at P in the direction v is 
defined as the sum of the sectional curvatures at P relative to the planes determined by v and the N — 1 directions 
orthogonal to v, i.e., if {ei, . . . , eN-i,v = cat} is an orthonormal basis of TpM and TTj is the plane spanned by v and 
ei. 
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AT-l 



Kn{P,v)=J2K{P,7ri). (A44) 



i=l 

The scalar curvature 7?. at P is the sum of the N Ricci curvatures at P, 

N 



7^(P)=^i^fl(P,ei) . (A45) 

i=l 

In terms of the components of the curvature tensor, such curvatures can be defined as follows (in the following formulae, 
we drop the dependence on P, because it is understood that the components are local quantities). We first define the 
Ricci tensor as the two-tensor whose components, Rij , are obtained by contracting the first and the third indices of 
the Riemann tensor, 

R,j = R% ; (A46) 

then 

Kr{v) = Rijv'v^ . (A47) 

The right hand side of Eq. (A47) is called "saturation" of Rij with v. The scalar curvature can be obtained as the 
trace of the Ricci tensor, 

Tl = R\. (A48) 

In the case of a constant curvature — or isotropic — manifold, the components of the Riemann curvature tensor have 
the remarkably simple form 

Rijki = K {gikQji - gugjk) , (A49) 
where K is the constant sectional curvature, so that the components of the Ricci tensor are 

Rij=Kgij, (A50) 
and all the above defined curvatures are constants, and are related by 



4. The Jacobi equation 

In this subsection we give a derivation of the Jacobi equation, already introduced in the main text as Eq. (27). We 

will proceed as follows: first, we will define the geodesic separation vector field J, then we will show that the field J 
is actually a Jacobi field, i.e., obeys the Jacobi equation. 

Let us define a geodesic congruence as a family of geodesies {Jt{s) = 7(s, r) | r e R} issuing from a neighborhood I 
of a point of a manifold, smoothly dependent on the parameter r, and let us fix a reference geodesic 7(5, tq). Denote 
then by 7(5) the vector field tangent to 7 in s, i.e., the velocity vector field whose components are 

y = f. (A52) 

and by J{s) the vector field tangent in tq to the curves 7s (r) for a fixed s, i.e., the vector field of components 

r = ^. (A53) 
ctr 

The field J will be referred to as the geodesic separation field, and measures the distance between nearby geodesies, 
as is shown in Fig. 28. 
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FIG. 28. Pictorial description of the definition of the geodesic separation vector J. From Ref. [22]. 

Let us now show that J is a Jacobi field. First of all, we notice that the field J commutes with 7, i.e., [7, J] = 0. In 
fact, from the definition of the commutator (Eq. A9) and from the definitions of J, Eq. (A53), and of 7, Eq. (A52), 
we have 

ds dx^ dr dx^ 

_ dr _ dY 

ds dr ' 



(A54) 



and using again Eqs. (A53) and (A52), we find that 



9J* d dx^ d dx^ 97* 
ds ds dr dr ds dr 



(A55) 



so that [7, J] = 0. Now, let us compute the second covariant derivative of the field J, V? J. First of all, let us recall 
that our covariant derivative comes from a Levi-Civita connection, which is symmetric (see Eq. (A28)), so that 

V^J-Vji=[i,J], (A56) 

and having just shown that [7, J] = 0, we can write 

J = Vj7 • (A57) 

Now, using this result, and the fact that V77 = because 7 is a geodesic, we can write 

= VjVjJ = V^Vj7 = [V^, Vj]7 , (A58) 

from which, using the definition of the curvature tensor (Eq. (A37)) and, again, the vanishing of the commutator 
[7, J], we get 

V? J = i?(7, J)7 , (A59) 
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which is nothing but the Jacobi equation (27), written in compact notation. 

It is worth noticing that the normal component J± of J, i.e., the component of J orthogonal to 7 along the geodesic 
7, is again a Jacobi field, since we can always write J = J± + A7: one immediately finds then that the velocity 7 
satisfies the Jacobi equation, so that J± must obey the same equation. This allows us to restrict ourselves to the 
study of the normal Jacobi fields, as we have already done in the main text. 



APPENDIX B: SUMMARY OF ELEMENTARY MORSE THEORY 



The purpose of this Appendix is to recall the main ideas and concepts of Morse theory which are relevant for the 
main text of the paper. For a more elaborate discussion we refer the reader to Refs. [118-120]. 

Morse theory, also referred to as critical point theory, links the topology of a given manifold M with the properties 
of the critical points of smooth (i.e., with infinitely many derivatives) functions defined on it. Morse theory links local 
properties (what happens at a particular point of a manifold) with global properties (the topology, i.e., the shape, of 
the manifold as a whole). Two manifolds M and M' are topologically equivalent if they can be smoothly deformed 
one into the other: a tea cup is topologically equivalent to a doughnut, but it is not topologically equivalent to a ball. 
In fact a ball has no holes, while both a tea cup and a doughnut have one hole. To define precisely what a "smooth 
deformation" is, one has to resort to the notion of a diffeomorphism. A diffeomorphism is a smooth one-to-one map, 
whose inverse is smooth. Then M can be smoothly deformed into M' if there exists a diffeomorphism tp which maps 
M into M' = ip{M). If such a diffeomorphosm exists, we say that M and M' are then diffeomorphic. Thus the notion 
of "topological equivalence" we referred to has now a precise meaning. 

For the sake of simplicity, we shall consider only compact, finite-dimensional manifolds: most of the results can 
be extended not only to noncompact manifolds, but also to infinite-dimensional manifolds modeled on Hilbert spaces 
(see [120]). 

The key ingredient of Morse theory is to look at the manifold M as decomposed in the level sets of a function /. 
Let us recall that the a-level set of a function / : M 1— > R is the set 

f-\a) = {xeM:f{x)=a}, (Bl) 

i.e., the set of all the points x G M such that f{x) = a. Now, M being compact, any function / has a minimum, /min, 
and a maximum, /max, so that 

/min < fix) < /ma. Vx € M . (B2) 

This means that the whole manifold M can be decomposed in the level sets of /: in fact, one can build M starting 
from /~^(/min) and then adding continuously to it all the other level surfaces up to /~^(/max)- To be more precise, 
one defines the "part of M below a" as 

Ma = r'{-oo,a] = {x€M : f{x)<a} , (B3) 

i.e., each Ma is the set of the points x <E M such that the function f{x) does not exceed a given value a; as a is varied 
between /min and /max, Ma describes the whole manifold M. 

For our purposes, we need to restrict the class of functions we are interested in to the class of Morse functions, 
which are defined as follows. Given a manifold M of dimension n and a smooth function / : M R, a point Xc (z M 
is called a critical point of / if df{xc) = 0, while the value f{xc) is called a critical value. The function / is called a 
Morse function on M if its critical points are nondegenerate, i.e., if the Hessian matrix of f at Xc, whose elements in 
local coordinates are 

= ' ^^^^ 

has rank n, i.e., has only nonzero eigenvalues. This means that there are no directions along which one could move 
the critical point, so that there are no lines (or surfaces, or hypersurfaces) made of critical points. As a consequence, 
one can prove that the critical points Xc of a Morse function, and also its critical values, are isolated. It can be proved 
also that Morse functions are generic: the space of the Morse functions is a dense subset of the space of the smooth 
functions from M to R. A level set f~^{a) of / is called a critical level if a is a critical value of /, i.e., if there is at 
least one critical point Xc S f~^(a). 

The main results of Morse theory are the following: 
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1. If the interval [a, 6] contains no critical values of /, then the topology of f~^[a,v] does not change for any 
V € {a,b]. This result^^ is sometimes called the non-critical neck theorem. The reason for this terminology will 
be made clear in the following. 

2. If the interval [a, b] contains critical values, the topology of f~^[a, v] changes in correspondence with the critical 
values themselves, in a way which is completely determined by the properties of the Hessian of / at the critical 

points. 

3. Some topological invariants of M, i.e., quantities that are the same for all the manifolds which have the same 
topology as M, so that they characterize unambigously the topology itself, can be estimated and sometimes 
computed exactly once all the critical points of / are known. 

Without giving explicit proofs, which can be found in Ref. [120] , let now us discuss in more detail items 1-3 above. 

1. The non-critical neck theorem 

If there are no critical values in the interval [a,b], there exists a diffeomorphism which sends f~^[a,b] into the 
Cartesian product f~^{a) x [a,b]. This means that the shape of f~^[a,b] is that of a multi-dimensional cylinder, or 
a neck (from which the name "non-critical neck"), if f~^(a) is simply connected, because the Cartesian product of 
a circle and an interval is a cylinder. This might be better understood with the aid of a two-dimensional example. 
Suppose that M is two-dimensional, and that the level set f~^{a) is topologically equivalent to a circle (see Fig. 29). 




FIG. 29. A non-critical neck. The lines with the arrows are the flow lines of V/, and the ellipses are the level sets of /. 

Then one can construct a diffeomorphism explicitly as the flow of the gradient vector field of /, V/, whose flow 
lines are orthogonal to the level surfaces of / and are depicted as the lines with the arrows in Fig. 29. This flow has 

no singularities if there are no critical values of /, so that the level set f~^{a) is transported up to f~^{b) along the 
flow lines of V/ without changing its topology. As a consequence, 

r'[a,b] « r\a) X [a,b] « r\b) X [a,b] , (B5) 

where "x ~ y" must be read as "x is difii'eomorphic to y" . 



We note that this result is valid even if / is not a Morse function; it is sufficient that it is a smooth function. 
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2. Critical points and topology changes 



In the neighborhood of a regular point P, N{P), there always exists a coordinate system such that / can be written 
as its first-order Taylor expansion^^, setting the origin of such coordinates in P, in the form 



fix) = m + ^x' + 



\/x e N{P) . 



(B6) 



Geometrically, this means that in the neighborhood of a regular point the level sets of / look like hyperplanes in R", 

because they arc the level sets of a linear function. 

But what if P is a critical point of /? A fundamental result by M. Morse, called the Morse lemma, is that if / is a 
Morse function then there always exists in N{P) a coordinate system (called a Morse chart) such that / is given by 
its second-order Taylor polynomial: 



/(^)=/(0) 



dx^dxi 



x'x^ + • 



Va; e N{P) . 



(B7) 



With a suitable rotation of the coordinate frame, {a;'} i— > {?/'}, the expansion (B7) can always be reduced to the 
canonical diagonal form 



f{y) = /(o) 



k 

i=l i=fe+l 



Vu e N{P) . 



(B8) 



Close to P, the level sets of / arc the level sets of a quadratic function, so that, geometrically, they arc non-dcgencratc 
quadrics, like hypcrboloids or ellipsoids, which become degenerate at P. The number of negative eigenvalues of the 
Hessian matrix, fc, is called the index of the critical point. Passing through the critical level, the shape of the level 
sets of / changes dramatically, in a way that is completely determined by the index k. Some examples in two and 
three dimensions are given in Fig. 30. 






4 / 




This follows from the implicit function theorem. 
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FIG. 30. Some examples of e-level sets near a critical point (the critical value of the function is set to 0). Upper left: n = 2, 
critical point of index k — 2; upper right: n = 2, critical point of index k = 1; lower left: n = 3, critical point of index k = 0; 
lower right: n = 3, critical point of index k — 2. 

The change undergone by the submanifolds Ma as a critical level is passed is described using the concept of 
"attaching handles". A /c-handle iJ^*^^ in n dimensions {0 < k < n) is a product of two disks, one fc-dimensional [D^) 
and the other (n — fc)-dimcnsional {D'^~^): 

ff(k) = £)fc X D"-*^ . (B9) 

In two dimensions, we can have either 0-handles, which are 2-dimensional disks or 1-handles, which are the product 
of two 1-dimensional disks, i.e., of two intervals, so that they are stripes, or 2-handles, which are again 2-dimensional 
disks (Fig. 31). In three dimensions, we have 0-handles which are solid spheres, 1-handles which are the product of a 
disk and an interval, so that they are solid cylinders, 2-handles which are the same as 1-handles, and 3-handles which 
are the same as 0-handles (Fig. 31). 




D(2) 1^0) jj(2) 

FIG. 31. Two-dimensional handles: ff*-"^ is the product of a 0-disk (a point) and a 2-disk, so that it is a 2-disk; Z/^^' is the 
product of two 1-disks, i.e., of two intervals, so that it is a strip; IS agam a 2-disk as IS. 
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FIG. 32. Three-dimensional handles: is the product of a 0-disk (a point) and a 3-disk (a ball), so that it is a ball; H^^' 

is the product of a 1-disk (an interval) and a 2-disk, so that it is a tube; ' is as and i^f^' IS as 

In more than three dimensions it is difficult to visualize handles: however, there is still the duality of the n ~ 2 and 
n = 3 cases, i.e., k and n — k handles are topologically equivalent. 

Having defined handles, we can state the main result of Morse theory as follows. 

Theorem. Suppose that there is one (and only one) critical value c in the interval [a,b], and that it corresponds to 
only one critical point of index k. Then the manifold Mi, arises from Ma by attaching a k -handle, and the transition 
occurs precisely at the critical level c. Everything goes in the same way if there are m > I critical points, with 
indices ki, . . . , on the critical level f^^{c); in this case Mi, arises from Ma by attaching m disjoint handles of types 
k\ , . . . , kjYi . 

Let us see how this works in a simple example. Consider as our manifold Af a two-dimensional torus standing on 
a plane (think of a tyre in a ready-to-roll position), and define a function / on it as the height of a point of M above 
the floor level. If the z-axis is vertical, / is the orthogonal projection of M onto the z-axis. Such a function has four 
critical points, and the corresponding four critical levels of /, which will be denoted as co,ci,C2,C3, respectively, are 
depicted in Fig. 33. We can build our torus in separate steps: each step will correspond to the crossing of a critical 
level of /. The steps are pictorially described in Fig. 34. As long as a < 0, the manifold Ma is empty. At a = cq = 
we cross the first critical value, corresponding to a critical point of index 0. This means that we have to attach a 
0-handle (a disk) to the empty set. Any Ma with < a < ci is diffeomorphic to a disk, as we can see by cutting a 
torus at any height between and ci and throwing away the upper part. At C2 we meet the second critical point, 
which now has index 1, so that we have to attach a 1-handle (a stripe) to the previous disk, obtaining a sort of a 
basket. Such a basket can be smoothly deformed into a U-shaped tube: in fact if we cut a torus at any height between 
ci and C2 and we throw away the upper part, we get a U-shaped tube. The third critical point C2 is again a point of 
index 1, so we have to glue another stripe to the tube. What we obtain can be smoothly deformed into a full torus 
with only the polar cap cut away from it. The last critical point has index 2, so that the crossing of it corresponds to 
the gluing of a 2-handle (a disk), which is just the polar cap we needed to complete the torus. 
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FIG. 33. The critical points and critical levels of the height function on a two-dimensional torus. 




FIG. 34. The building of a two-dimensional torus by attaching handles, (a) Attaching a 0-handle. (6) attaching a 1-handle. 
(c) Attaching a 1-handle again, (d) attaching a 2-handle to complete the torus. The symbol "~" means "is diffeomorphic to" . 

3. Topological invariants 

Morse theory can be used also to give estimates, and sometimes to compute exactly, some topological invariants 
of our manifold M . For a (two-dimensional) surface, a very important topological invariant is the genus g, which 
equals the number of holes of the surface. The generalization to n dimensions of the genus is given by the set of the 
Beta numbers bk{M), < k < n, which are the fundamental topological invariants of an n-dimensional manifold, and 
completely describe its topology. For hyperspheres, all the Betti numbers are zero. Using the Betti numbers we can 
obtain another topological invariant, the Euler characteristic x{M), which is nothing but the alternating sum of the 

n 

X(M) ^ Y.{-l)%^{M) . (BIO) 

k=l 

For (two-dimensional) surfaces, x = 2 — 2(7 holds. 

Now let us consider a generic Morse function f on M and let us define the Morse numbers ^k{M) as follows: is 
the total number of critical points of / on M which have index k. It turns out that the Morse numbers of a manifold 
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are upper bounds of the Betti numbers, i.e., the following (weak) Morse inequalities hold: 



hk{M)<txk{M) 



k = 0, . . . ,n . 



(Bll) 



Actually, a result stronger than Eqs. (Bll) holds, which states that alternate sums of two, three, four, ... subsequent 
Bctti numbers arc bounded from above by the alternate sums of the corresponding Morse numbers. Starting from 
this result one can prove the following identity 



and this provides a way of computing exactly the Euler characterstic of a manifold once all the critical points of a 

Morse function are known. 

Among all the Morse functions on a manifold M, there is a special class (called perfect Morse functions) for which 
the Morse inequalities (Bll) hold as equalities. Perfect Morse functions characterize completely the topology of a 
manifold. It is possible to prove that the height function on the torus we considered above is a perfect Morse function 
[120]. However, there are no simple general recipes to construct perfect Morse functions (this is actually an active 
area of research). 



For a long time the equations of Newtonian mechanics have been the paradigm of classical determinism. Only quite 
recently has it been realized that "determinism" and "predictability" are far from being the same concept, and that 
predictability also requires the stability of the solutions of the dynamical differential equations. Determinism implies 
that, once an initial condition is given, the trajectory is uniquely determined for all the forthcoming times; stability 
means that two initially close trajectories will remain close in the future (more precisely, their distance will grow 
slower than a power of the time). If this is not true, it becomes impossible to predict the evolution of a system even 
for very small times, as explained, for instance, by R. P. Feynman [121]: 

It is true classically that if we knew the position and the velocity of every particle in the world, or in a 
box of gas, we could predict exactly what would happen. And therefore the classical world is deterministic. 
Suppose, however, that we have a finite accuracy and do not know exactly where just one atom is, say to 
one part in a billion. Then as it goes along it hits another atom, and because we did not know the position 
better than one part in a billion, we find an even larger error in the position after the collision. And that 
is amplified, of course, in the next collision, so that if we start with only a tiny error it rapidly magnifies 
to a very great uncertainty. To give an example: if water falls over a dam, it splashes. If we stand nearby, 
every now and then a drop will land on our nose. This appears to be completely random, yet such a 
behavior would be predicted by purely classical laws. The exact position of all the drops depends upon 
the precise wigglings of the water befor it goes over the dam. How? the tiniest irregularities are magnified 
in falling, so that we get complete randomness. Obviously, we cannot predict the position of the drops 
unless we know the motion of the water absolutely exactly. 

Speaking more precisely, given an arbitrary accuracy, no matter how precise, one can find a time long 
enough that we cannot make predictions valid for that long a time. Now the pint is that this length of 
time is not very large. It is not that the time is millions of years if the accuracy is one part in a billion. 
The time gose, in fact, only logarithimically with the error, and it turns out that in only a very, very 
tiny time we lose all our information. If the accuracy is taken to be one part in billions and billions and 
billions no matter how many billions we wish, provided we do stop somewhere then we can find a 
time less than the time it took to state the accuracy — after which we can no longer predict what is going 
to happen! 

As long as nonlinear dynamical systems are considered, stability is the exception rather than the rule. Even if this 

relies - - at least from a conceptual point of view — upon mathematical results which have been known since the end 
of the last century, its importance has only been completely realized with the aid of a new and powerful approach: 
numerical simulation. The very complicated structure of some trajectories which can arise in nonlinear dynamical 
systems was discovered by Poincare [1] in the late XIX century, but the physics community became fully aware of 
the existence and of the meaning of these structures only they were visualized by computer simulation in the work of 
Henon and Heiles [59]. 



n 




(B12) 



fe=i 



APPENDIX C: CHAOS IN HAMILTONIAN DYNAMICAL SYSTEMS 
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The instability wc are referring to is known as intrinsic stochasticity of the dynamics, or "deterministic chaos" . 
These terms mean that the dynamics, being completely deterministic, yet exhibits some features that make it in- 
distinguishable from a random process. The characteristic feature of a chaotic system, which is at the basis of the 
unpredictability of its dynamics, is the sensitive (exponential) dependence on initial conditions: the distance between 
two trajectories which originate in very close-by points in phase space grows exponentially in time so that the system 
looses the memory of its initial conditions. Regular dynamics, i.e., quasiperiodic motion, is — as far as conservative 
systems are considered — a "weak" property, because it is destroyed by very small perturbations of the system. On 
the contrary, chaos is a strong property, because given a dynamical system where chaos is present, in many cases it 
will be present even after the system has been subjected to significant perturbations [122]. 

Here we recall briefly the main concepts of the theory of Hamiltonian dynamical systems which are necessary for the 
understanding of the material on chaos presented in this Report. The main goal of this Appendix is then to provide 
the reader with a definition of the Lyapunov exponents and of a motivation for the introduction of these quantities 
as a "measure" of chaos in a dynamical system. 

A very good introduction to the subject is given in Lichtenberg and Lieberman's classic book [122], and, at a more 
pedagogical level, in Tabor's [123] and Ott's textbooks [124]. An interesting selection of reprints can be found in 
MacKay and Meiss [125]. We assume the reader is familiar, at least at a basic level, with the concepts of ergodicity 
and mixing. A discussion on these topic would be far beyond the scope of the present Report; a good introduction 
on these topics can be found in any of the references just mentioned above. 



Throughout the paper wc have been concerned with Hamiltonian dynamical systems with a large number of degrees 
of freedom. However, the main features of chaos can be better appreciated starting with an example of a system with 
only one degree of freedom, subjected to an external perturbation: the forced pendulum. Although the behavior of 
many-degrec-of-frccdom is much more complicated, nevertheless some of the essential features of chaos are already 
present in this simple example. 

The forced pendulum is a system obeying the follwing equation of motion: 



The phase space of the forced pendulum is three-dimensional, because in addition to the coordinate q and to the 
momentum p = q, one has to take into account also the time t, because the forcing term on the right-hand side of Eq. 

(CI) depends explicitly on time. The forced pendulum is, however, a Hamiltonian system, and since the dimension of 
the phase space of autonomous Hamiltonian dynamical systems is 2iV, where N is the number of degrees of freedom, 
it is customary to refer to systems like the forced pendulum as to systems with "1.5 degrees of freedom". 
As long as e = 0, the system obeying Eq. (CI) is a simple pendulum, and its Hamiltonian 



is an integral of motion, so that its value, the energy E, is a constant of the motion and the system is integrable, as every 
one-degree-of- freedom autonomous Hamiltonian systems is. The word "integrable" is used here in a wider sense than 
its immediate meaning "such that the equations of motion can be solved" ; a Hamiltonian system is integrable when 
it has a sufficiently large number of integrals of motion {N, for an autonomous system with N degrees of freedom), 
such that its trajectories do not explore the whole phase space, but are confined to lower-dimensional subsets called 
invariant tori^^ . When N = 1, each invariant torus coincides with a trajectory. Some of these are depicted in Fig. 
35. We remark that there arc two distinct classes of trajectories: oscillations, which correspond to bounded motions, 
and rotations, which are unbounded, the two classes are separated by a curve called the separatrix. The separatrix 
is the trajectory pursued by the pendulum when it starts precisely at the unstable equilibrium point (p, q) = (0, ±7r) 
with E = 1, i.e., just the energy that is required to come back to the same point (note that q = n and q = —n must 
be identified). The motion on the separatrix requires an infinite amount of time. 



The origin of the term is as follows: such subsets are invariant because if a trajectory starts on one of them, it remains 
there forever; they are called tori, because they are topologically equivalent to multi-dimensional tori. 



1. A simple example of chaotic dynamics: the perturbed pendulum 



^' -|- sin g = e sin(g — ivt) . 



(CI) 




cosq 



(C2) 
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FIG. 35. Phase-space trajectories of a simple pendulum. The closed curves are the oscillations, the curves above and below 
the separatrix S are clockwise and counterclockwise rotations, respectively. 

But what happens if e ^ 0? Once e ^ 0, no matter how small, the system (CI) is no longer integrable, and the 

separatrix, which was a unique curve in the e = case, splits into two distinct invariant curves. These curves must 
intersect transversally each other infinitely many times, as Poincare showed for the first time [1]. These intersections 
are referred to as homoclinic intersections, and force the trajectories to fold themselves giving rise to a very complicated 
structure: in Poincare's own words [1] , 

these intersections form a sort of texture, or of a net whose meshes are infinitely tight; each of these 
two curves can never intersect itself, but has to fold in a complicated way as to intersect all the meshes of 
the net an infinite number of times. One is amazed by the complexity of this picture, which I do not even 
attempt to draw. 

As a consequence of the presence of these intersections, in a neighborhood of the region in phase space which was 
occupied by the separatrix in the integrable case, a so-called chaotic sea suddenly appears. The chaotic sea is the 
region irregularly filled by dots in Fig. 36, where a two-dimensional section^^ of the 3-d phase space of the system 
is shown. If we now follow the evolution of two intially close points in the chaotic sea, we see that their separation 
grows exponentially in time, so that the dynamics in the chaotic sea is unpredictable. 



^'^This section has been obtained as a stroboscopic Poincare section [123], so that each point on the plot corresponds to an 
intersection of a trajectory of the system with the planes t = 2n7r/w. 
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FIG. 36. Section of the phase space of a perturbed pendulum, showing the appearance of chaotic seas close to the separatrix 
of the unperturbed system (the solid hnc in the figure). The dots are obtained from a single trajectory issuing from a point 
very close to the unperturbed separatrix. The amplitude of the perturbation is e = 10"*. 

The appearence, in phase space, of irregular regions like the chaotic sea could justify by itself the use of the term 
"chaotic dynamics". However, there are also other properties of the dynamics described by Eq. (CI) which justify 
the use of such a term. For example, if we introduce a symbolic coding of the dynamics in which the symbol is 
associated with each passage through the point g = with g > and the symbol 1 to each passage through the same 
point with g < 0, then given any bi-infinite sequence of zeros and ones, for example generated by coin tosses, this 
sequence corresponds to a real trajectory of the system (CI). The motion of the system, though deterministic, is thus 
indistinguishable from a random process. 

We can intuitively understand the origin of such a behaviour if we think that when the phase point is on a trajectory 
very close to the separatrix, an infinitesimal variation may qualitcdively change the character of the motion (e.g.. from 
oscillations to rotations). This is an example of the sensitive dependence on initial conditions that, in general, leads 
to the exponential separation of initially close-by trajectories. 

This example is extremely simple but contains the essential features of the problem. In fact even the case of a 
Hamiltonian system with N degrees of freedom can be treated in a similar way and shows analogous, though much 
more complicated, results. In that case, there is no need of an external forcing to get chaos, for as A/' > 1 an 
autonomous, nonlinear Hamiltonian system is generically nonintegrable (the integrable systems being a very small 
subset of all the possible systems). However, even in the simple low-dimensional cases, by means of concepts like 
homoclinic intersections, it is possible only to give a qualititative description of the onset of chaos, but a quantitative 
description of the stochastic regions is impossible, i.e., it is impossible to compute how fast two initially close-by 
points will separate. For A/'-dimensional systems, the situation is obviously even worse: there exists a method, which 
is a generalization of Poincare's method, obtained by Mel'nikov [126] and later by Arnol'd [127], which allows one to 
show the existence of homoclinic intersections near separatrices for very small perturbations even for large systems, 
but again there is no possibility of describing quantitatively the stochastic regions. 

To obtain quantitative informations on chaotic dynamics we must introduce the concept of Lyapunov exponent. 

2. Lyapunov exponents 

We now give a definition and an explanation of the Lyapunov exponents. Our discussion will be aimed at showing 
how to define and compute the Lyapunov exponents for a dynamical system which is defined by a system of ordinary 
differential equations, i.e., for a flow, because Hamiltonian dynamical systems belong to this class. For a more general 
discussion of Lyapunov exponents, see Eckmann and Ruelle [128]. 

From a physical point of view, given a trajectory of a dynamical system, it is important to find answers to the 
following questions: Is the trajectory chaotic? And if so, how strong chaos is, or how fast do two initially close-by 
points separate in phase space, i.e., how long should we wait until the system exhibits its chaotic feature? The 
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concept of Lyapunov exponents is introduced to answer these questions, since Lyapunov exponents are defined in 
order to provide an average measure of the rate of exponential divergence of nearby orbits in phase space, which is 
the distinctive feature of chaos. 

Lyapunov exponents are defined for a given trajectory of a dynamical system, and this allows us to give a definition 
of a chaotic trajectory as follows: a trajectory is said to be chaotic if its (largest) Lyapunov exponent is positive. 

We now give a definition of the Lyapunov exponents. Let us consider a dynamical system whose trajectories in an 
n-dimensional phase space M are the solutions of the following system of ordinary differential equations: 

Xi = Xi{xi, ...,Xn) , 

: : (C3) 

Xn — -^71 {xi , . . . , Xn ) , 

If we denote by x{t) = {xi{t), . . . ,Xn{t)) a given trajectory whose initial condition is x(0), and by y{t) another 
trajectory which is initially close to x{t), and we denote by ^{t) the vector 

m = Vit) - --rit) , (C4) 

then the evolution of ^ describes the separation of the two trajectories in phase space. The vector ^ is assumed to 
obey the linearized equations of motion, because it is assumed to be initially small. These equations are, as can be 
shown by inserting Eq. (C4) into the equations of motion (C3) and expanding in a power series up to the linear terms. 



1 V 9Xj 



x{t) 



(C5) 



dxi 



j=l ^ -J / x{t) 



and are referred to as the tangent dynamics equations^^ , which we already wrote in the main text in the particular 
case of a standard Hamiltonian system (see Eq. (25)). Note that (C5) is a system of linear differential equations, 
whose coefficients, however, depend on time. According to the definition (C4), the norm |^| of the vector ^, i.e.. 



\m\ = 



1/2 



(C6) 



measures the distance of the two trajectories as a function of t. If the trajectory x{t) is unstable, all its perturbations 
grow exponentially, so that oc exp(At). If the elements of the Jacobian matrix dXi/dxj, which are the coefficients 
of the linear equations (C5), were either constant or periodic, it would be possible to solve the system, but, since the 
Jacobian matrix depends on the trajectory a;(f), its entries are in general neither constant nor periodic, so that the 
rate of exponential divergence varies with time. Therefore, one introduces an asymptotic rate of exponential growth 
of ^ as the Lyapunov exponent A 

A = hm - log — — , (C7) 

t^oo t |C(0)| 

which measures the degree of instability of a trajectory: if A is positive, the trajectory is unstable with a characteristic 
time A~-^. In principle A depends on both the initial values of x, x(0), and of ^, ^(0). However, Oseledec [129] has 
shown that the limit (C7) exists, is finite and can assume only one of the n values 

Ai < As < • • • < A„ . (C8) 

The set {Aj} is the called Lyapunov spectrum . The exponent A defined in (C7) takes the n diff'erent values of the 
spectrum as the initial condition ^ in the tangent space Tx{o)M is varied; the latter admits a decomposition in linear 

subspaces. 



^This notation follows from that the dynamics of the vector ^ takes place in the tangent space T^^t^M of the phase space M. 
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T^(o)M = ^1 e £^2 e • • • e , (C9) 

and each is associated with the corresponding subspace Ei, in that a vector ^(0) G Ei will exponentially grow with 
the exponent Aj. If there exists on the phase space M a probability measure //, which is ergodic and invariant for 
the dynamics on M, then the numbers Aj do not depend on the initial condition x(0), apart from a possible subset of 
initial conditions of measure zero with respect to /i. 

In practice, the evolution of the norm of a tangent vector is sensitive only to the first — the largest — exponent, 
because a generic initial vector ^(0) will have a nonvanishing component in the E\ subspace, so that the largest 
exponent Ai will always dominate in the long-time limit: choosing ^ at random with respect to a uniform distribution 
we have A = Ai with probability one. This means that Eq. (C7) provides a practical definition for the largest Lyapunov 
exponent Ai, which we have always denoted simply by A in the main text. 

Let us now apply the above to a standard Hamiltonian system, whose Hamiltonian is of the form (1); the dimension 
of the phase space is n = 2iV, and the equations of motion (C3) are now Hamilton's equations 

qi = Pi , 

qN = Pn , 
dV 

Pi = (CIO) 
oqi 

dV 

Pn = -- 



dqN 

and also the linearized dynamics (C5) can be cast in the canonical form 

6 = Cjv+1 , 



N 



^2^ = -Ef7Sr)..^^- 



KdqNdqjy^^^^ 

This equation was already introduced as Eq. (25) in Sec. Ill, and is usually referred to as the tangent dynamics 
equation for Hamiltonian systems. To measure the largest Lyapunov exponent A in a numerical simulation, one 
integrates numerically both Eqs. (CIO) and (Cll), and then makes use of the definition (C7), which can be rewritten, 
in this case, as 



A = lim - log 



eiit)+---+eNit)+ei{t)+---+eNit) 



?f(o) + --- + e^(o) + ?f(o) + --- + S^(o) 



1/2 

— - (C12) 



where we have used that = ^i+N (see Eq. (Cll)). More precisely, in a numerical simulation one uses the discretized 
version of Eq. (C12), i.e., 

A = hm — > — log ' ,, , (C13) 

m^oom'^At ^ |C(«Ai)| ^ ' 

where, after a given number of time steps At, the value of |^| has to be renormalized to a fixed value, in order to 
avoid overflow [130]. 

The definition (C7) does not allow one to measure the other exponents of the Lyapunov spectrum. To measure 
them, one has to observe that they can be related to the growth of volumes in the tangent space. A two-dimensional 
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area V2 in the tangent space, spanned by two linearly independent tangent vectors ^^^^ and ^^^^ will expand according 
to 

FaW ocexp[(Ai +A2)t] , (C14) 

a three-dimensional volume, as 

Vsit) oc exp[(Ai + A2 + Xsjt] , (C15) 
and so on, so that, choosing k <n linearly independent and normalized vectors ■ ■■ , ^^^'^ S T^M we obtain 

k 

lim - \og\e\t)he''\t)A---A =^A, , (C16) 

i=l 

Therefore the algorithm (C13) can be generalized to obtain an algorithm to compute numerically the whole Lyapunov 
spectrum [130]. However, such a computation is very hard when the number N is large. 

The sum of all the n Lyapunov exponents in the Lyapunov spectrum, X^"_iAj, measures the expansion rate of 
n-volumes in phase space. Therefore, for a Hamiltonian system, 

2JV 

^Ai = 0, (C17) 

i=l 

because volumes in phase space are conserved. In addition, for Hamiltonian systems a result stronger than (C17) 
holds, i.e., there is a symmetry in the Lyapunov spectrum such that 

A, = — A2jv-j+i • 

Equation (C18) for Hamiltonian systems is a consequence of the symplectic structure of the Hamilton's equations 
[131], however it has been recently generalized to a class of non-Hamiltonian systems [132]. 

The numerical integration of the Eqs. (C5) and the consequent measure of A — or of the spectrum {Ai} when 
it is possible in practice is the standard technique to characterize Hamiltonian chaotic dynamics. An operative 
definition of a chaotic dynamical system can be stated as follows: a system is chaotic if it has at least one positive and 
one negative Lyapunov exponent. In fact this ensures that the system shows (almost everywhere with respect to the 
ergodic measure /i used to define the Lyapunov exponents) the distinctive features of chaos as described in the example 
of the forced pendulum. In fact, the presence of a positive exponent ensures the presence of a exponential divergence 
of nearby orbits, and the presence of a negative one ensures that they also fold and mix in a very complicated way, so 
that they can produce those structures we refrred to as "chaotic seas" . However, as long as autonomous Hamiltonian 
systems are considered, the anti-symmetry of the spectrum (C18) ensures that the presence of a positive exponent 
implies the presence of a negative one with the same absolute value, so that a single (the largest) positive exponent is 
sufficient to have chaos; on the contrary, if the largest exponent vanishes the dynamics will be regular. These facts, 
together with that the largest Lyapunov exponent A measures the smallest instability time scale, show how natural 
the use of the value of A is to measure chaos in such systems. 

It is important to specify with resepct to what invariant ergodic measure ji the Lyapunov exponents are defined: 
this may be also a J-measure concentrated on a single trajectory, in which case we could speak of a chaotic trajectory 
rather than of a chaotic system. In Hamiltonian systems with a large number of degrees of freedom we expect the 
microcanonical measure of the chaotic regions to be overwhelmingly larger than the measure of the regular regions; 
the existence of these regular regions is ensured — at least as long as the system is not too far from an integrable limit 
— by the Kol'mogorov-Arnol'd-Moser (KAM) theorem [56]. However, from a practical point of view the measure 
relevant for the definition of the Lyapunov exponent is indeed the microcanonical one. Numerical experiments are 
in agreement with this expectation for large systems, since no relevant dependence of the Lyapunov exponent on the 
initial conditions has been detected, and this is the reason why in the main text we have never referred explicitly to 
any dependence of A on /x, treating the Lyapunov exponent as any other "thermodynamic" observable. Nevertheless 
for small systems (especially N ~ 2 which is the best known case) the simulations show that the measure of the 
chaotic regions may be very small in a very large energy range, so that in that case one cannot speak of a truly chaotic 
system but only of a system in which chaotic and regular regions are simultaneously present (these systems are often 
referred to as mixed systems, as they are in between completely chaotic and regular ones). 

Since we are interested in large systems, up to the thermodynamic limit, a number of questions naturally arises: 
what is the behavior of the Lyapunov exponents as n increases; does a thermodynamic limit exist for the Lyapunov 
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spectrum, etc. Numerical results [133] have shown that as n — > oo the Lyapunov spectrum {Aj} appears indeed to 
converge to a well-behaved function 

X{x) = lim . (C19) 

n— ^oo 

The function A(a;) is a non-increasing function of a; € [0,1]. Some rigorous work in this respect has been recently 
done by Sinai [134] . The existence of a limiting Lyapunov spectrum in the thermodynamic limit has many important 
consequences that we will not review here; a good discussion can be found in Ref. [131]. We only want to remark here 
that the existence of a thermodynamic limit for the Lyapunov spectrum implies that the largest Lyapunov exponent 
is expected to behave as an intensive quantity as increases. 



APPENDIX D: THE STOCHASTIC OSCILLATOR EQUATION 



In the following we will briefly describe how to cope with the stochastic oscillator problem which we encountered 
in §IVA3. The discussion closely follows Van Kampen (Ref. [65]) where all the details can be found. 
A stochastic differential equation can be put in the general form 



F{x,x,x, . . . ,0) = 0, 



(Dl) 



where F is an assigned function and the variable is a random process, defined by a mean, a standard deviation and 
an autocorrelation function. A function ^(O) is a solution of this equation, if VO F(^(f2),f2) = 0. If equation (Dl) is 
linear of order n, it is written as 



u = A{t, 0)u 



(D2) 



where 



u = 



Us 



X \ 

X 
X 



(D3) 



and A is an n X n matrix whose elements A^,^{t) depend randomly on time. 

For the purposes of our work, we are interested in the evolution of the quantities u^u^, rather than of the u^'s 
themselves. The u^i/'s obey the differential equation 



dt 



{t){UkUx) ,] 

fe,A 



where 



(D4) 



(D5) 



However, both Eq. (D2) and (D4) have exactly the same form and can be solved using the same procedure, so that 
we will first illustrate such a procedure in general. Therefore in the following formulae, u refers to a vector whose 
components are either the u^'s or the u^^'s, and A denotes either the matrix A in Eq. (D2) or the matrix A whose 
elements are given by Eq. (D5), respectively. Then, we will apply this procedure to the case of the stochastic harmonic 
oscillator. 

Now, solving a linear stochastic differential equation means determining the evolution of the average of u{t), (u(i)), 
where the average is carried over all the realizations of the process. Let us consider the matrix A as the sum 



A{t, O) = Ao(i) + aAi(t, fl) 



(D6) 



where the first term is il-independent and the second one is randomly fluctuating with zero mean. Let us also assume 

that Ao is time-independent. If the parameter a. which determines the fluctuation amplitude is small we can 
treat Eq. (D2) by means of a perturbation expansion. It is convenient to use the interaction representation, so that 
we put 
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u(i) = exp(Aot)v(t) (D7) 

and 

Ai{t) = exp(Aot)v(f) exp(-Aot) . (D8) 

Formally one is then led to a Dyson expansion for the solution v(i). Then, going back to the previous variables and 
averaging, the second order approximation gives 

d f~^°° 

-{u{t)) = {Ao + a^ (Ai(t)exp(AoT)Ai(f-T))exp(-Aor)dT}(u(i)) . (D9) 

Let us remark that, if the stochastic process 11 is Gaussian, Eq. (D9) is more than a second order approximation: 
it is exact. In fact, the Dyson series can be written in compact form as 



{u{t))=T 



exp (^J^ A{t')dt' 



(u(0)) , (DIO) 



where T[- ■ ■] stands for a time-ordered product. According to Wick's procedure we can rewrite Eq. (DIO) as a cumulant 
expansion, and when the cuimilants of higher than the second order vanish (as in the case of a Gaussian process) one 
can easily show that Eq. (D9) is exact. 

We now apply this general approach to the case of interest for the main text, i.e., to the stochastic harmonic 
oscillator equation, which is the the second-order linear stochastic differential equation given by 

x + n(t)x = 0, (Dll) 

where il{t) is the random squared frequency, O = Oq + crj^r]{t), where fio is the mean of ^{t), an is the amplitude of 
the fluctuations, and r]{t) is a stochastic process with zero mean. In this case, Eq. (D2) has the form 



d / x\ ( 1 



dt\x ) ^ )\x 



X 



(D12) 



In particular, we are interested in obtaining the averaged equation of motion for the second moments. Using Eqs. 
(D5) and (D12), one finds that Eq. (D4) becomes: 





(D13) 



XX 





Like in Eq. (D6), the matrix A splits into: 

K{t) = Ao + anr?(f)Ai = | -2^0 | + CTor/(t) | -2 | , (D14) 

so that the equation for the averages becomes 

- {x?) ={Ao + 4 / (^(t)r?(i-t'))B(t')rft'} ki?) , (D15) 
^* \ (xx) J Jo y {xx) J 

where B(t) = Ai exp(Ao<)Ai exp(— Agt). 

When the process 'r]{t) is Gaussian and (5-correlated, Eq. (D15) is exact, and the integral can be computed explicitly: 
writing {r]{t)r]{t — t')) = T5{t'), where r is the correlation time scale of the random process, we obtain 



{Ao + ^B(0)}| (P) I . (D16) 



{xx) 

Prom the definition of B(f) it follows then that B(0) = Af , and by an easy calculations we find 




78 



Ao + oIt Al = 

which is the result used in 81 V A 3. 



2 \ 

air -2^0 (D17) 
-9.0 10/ 
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